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Numerical Calculation of Certain Definite 
Integrals by Poisson’s Summation Formula 


It is generally agreed that of all quadrature formulae, the trapezoidal rule, 
while being the simplest, is also the least accurate. There is, however, a rather 
general class of integrals for which the trapezoidal rule can be shown to be a 
highly accurate means of obtaining numerical values. Specifically, if the integrand 
is a periodic, even function with all derivatives continuous, then the integral over 
a period is given by the trapezoidal rule with high precision. Similarly, integrals 


of the form f ‘ e~* f(x)dx, where f(x) is an even, entire function of x can be 


0 
evaluated by this method. 
The basis for the above statements is the Poisson summation formula [6] 
which reads 


[soa =h {sL710) + f(a)) + = f(kh)} — 2 = (7) 
(1) zm i 


h=a/n, g(x) = ih f(t) cos (xt)dt. 


This formula may be interpreted as the trapezoidal rule with a correction. The 
magnitude of the correction will depend on the manner in which the Fourier 
coefficients, g(x), decrease with increasing x. Now if f(#) is periodic, of period a 
with all derivatives continuous, then it is known from principles of Fourier series 
that the coefficients in the series for f(#) tend to zero at least as fast as e~* for x 
sufficiently large. If this is the case, the first term in the correction to the trape- 


2x 
zoidal rule will be of the order of e *, which is small even when “hk” is compara- 
tively large. Thus, if it is possible to estimate in advance the order of magnitude 
of g(x) for large values of x, an interval of integration may be determined to give 
any desired degree of accuracy. 

If, in Poisson’s formula, a and m are made to tend to infinity in such a way that 
h remains fixed, there is obtained the result 


J seae = 2 {470 +E so} -2E (4) 


(2) y 
g(x) = f f() cos (xt)dt. 


Thus, if f(x) is an even function such that its Fourier transform g(x) is negligibly 
small for large x, the trapezoidal rule may again be used with an arbitrarily 
small error. 

We consider now some functions to which the above may be applied. 
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(1) The Bessel Function J,(2), of Integral Order. This function may be defined 
by the definite integral as follows: 


(3) Jp(s) = 2 [ cue cos (pit, 
where “‘p’’ is a positive integer. 
Thus 
f(t) = e* =" * cos (pt); 
(4) g(x) = {" e* 8 t cos (pt) cos (xt)dt 
0 


Il 


© Lie? Serp(2) + i” Jeg(s)], 


(x an integer). 




















To obtain an error estimate from (4), it is noted that if ‘‘x” is chosen suffi- 
ciently large, so that J.,(z) < Jz_,(z), then the order of magnitude of g(x) is 
essentially that of J,,(z), and this in turn may be made arbitrarily small by the 
proper choice of ‘‘x,”’ assuming “‘p” and “‘z’”’ fixed. Therefore, the order of magni- 
tude of the error, for which only the first term in the infinite series in (2) need be 
considered, is given by J,_,(z), with x = 2x/h. But by definition, h = x/m, and 
so the error estimate in terms of the number of subdivisions, », of the interval 
of integration is given by Jen_,(z). 

The value of “‘h’’ thus determined permits J,(s) to be approximated by a finite 
sum of sines and cosines for any value of z and p, and to any accuracy which may 
be required. We thus obtain from the Poisson summation formula a simple 
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analytic approximation to the Bessel function of the first kind of integral order, 
involving only elementary functions. Such a representation provides a means of 
generating these functions on an automatic digital computer using only subrou- 
tines for the trigonometric functions. In figure 1, the value of logis J:,(z) is 
plotted as a function of z for those values of m for which J;,(z) lies between 10-* 
and 10~'*. The use of this plot is as follows: Assume that it is desired to approxi- 
mate Jo(z) in the range 0<z<2, with a prescribed error of the order of 10-*. From 
the figure we see that when 2m = 12 (whence h = x/6) the error is of the desired 
order of magnitude. Therefore we find that, for the range 0 < z < 2, 


xJo(z) one [te* + et cos +/6 + e%# cos #/3 + eit cos x/2 


+ he-* + eit cos 5/6 + gt ene tose), 
or 


1 
(5) Jos) ~4| cos (2) + 2.008 (2) + 2.008 (42) +1], 
correct to at least eight decimal places. To obtain ten decimal place accuracy we 
need only increase the number of subdivisions from 6 to 8. For z = 1, formula 
(5) gives 

Jo(1) = .76519 76867 


which is correct to nine places with an error of one unit in the tenth. The values 
obtained by applying respectively Simpson’s rule and Weddle’s rule, with the 
same interval, are 

Jo(1)], = .76521 


Jo(1) he = .76572. 


With A = 2/15, Jo(z) is given correctly to ten places for all values of z up to 11 by 
6) J 2 (2) +20 | +os % | + 2c0s| 2 =| 
- me~oy " 15 8 15 


v an 
+ 2.cos| # cos | + 2.cos| scos | 


2 
+ 2.08 | s 08 | + 2008 008 7 | + 2.08 | 0082 | . 
3 5 15 
(2) The Modified Bessel Function I,(2). This function can be computed in a 
similar manner by replacing z by iz. J,(z) does not tend to zero quite as fast as 


does J,(z) when n — © but the difference is not appreciable if z is not large [12]. 
The approximation (5) gives 


(7) In(s) = 1 cosh (s) + 2 cosh (%:) + 2 cosh (3s) + | 
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to within 10-'° if ¢ < 1.2 and to within 10~’ if z < 2. For z = 1 we obtain 


From (7): Io(1) = 1.266065878. 
From Simpson’s Rule: Jo(1) = 1.266051. 
From Weddle’s Rule: Jo(1) = 1.26658. 


The correct value from [8] is 1.2660658778. 

The approximations obtained in this way may also be used to evaluate in- 
tegrals involving Bessel Functions. For example the functions considered by 
ScHwartz [10] defined as 


Je(h, x) = f * Jo(At) cos t dt, 
0 

Ji(h, x) = f * Jo(dt) sin t dt, 
0 


may be evaluated in this way. If x < 2, the approximation (5) gives, after carry- 
ing out the integrations, 





200s ( r«) 
cos (Nt) yess & 2 cos m 
Ce, 1— 2? 1 —})? 


8 sin (“3 ae) 
- see | (Ax) 4 2 = || 


1—» 1 — 32 1 — })? 


(8) J.(A, x) ~ ; 4 x | 








A similar approximation may be obtained for J,(A, x). 
(3) The Modified Functions of the Second Kind, K,(z). These are defined by 


(9) K,(s) = f * g-teoth s cosh (nx)de. 


0 


Hence 


g(x) = f 7 e~* 8h § cosh (mt) cos (xt)dt 
0 


0 
(10) = $[Ko+ie(s) + Knsic(2)]. 
Now 
2 — Je(#) — Ip(2) 
(11) 5 ee. Se 


Thus, using the approximation for |p|>> z: 


(2/2)? 
T'(p + 1)’ 


I,(%) ~ 
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(2)| ei In (2/2) giz In (2/2) 
TAR r(iit+n—ix) T(i+n+ ix) 


Hence if 
(i —m+ ix) = |T'(1 — 2 + ix)|e™, 


Tai ++ ix) = |T(1 +2 + ix)|et™, 
then, 


1 
(13) =| Knsie(2) + K_n+i2(2) | 


‘ (on sin (» —xIn (2)) 


~ sinh (xm) |r(1 — n + ix)| 


(sy nbec en) 


|r(1 +” + ix)| 


GG) 











( 








1 
pate Se” |T(1 + + ix)| 


But 


|P(1 + m + ix)|? = [n? + 27] (m — 1)? + 2°) --- 11 + JL? Tes) |? 


| T (ax) |? 





a ix) |2 = | 
(14) |r(i — + ix)| [@—i?+e)- +2] n> 
= |T'(ix)|?, n=1; 
== x*°|I' (sx) |*. n= 0. 
Therefore, since |T'(ix)|* = eee, 
x sinh (xx) 


. (5) 


Vx sinh (wx) | V[1 + x*] --- [n? + x7] 





|Knsie(2) + K—n+ie(2)|~ 











, 20+ 9642) Te Fe) 
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@) 4} 


form > 1; 





(15) | Kizie(2) + K-14ie(2)|~ 





+- 


= | 
xV¥1 + x? z 


Vsinh (xx) 


> 
Vx sinh (xx) 


The last result agrees with a more accurate expression for K;,(z) when x > z 
(cf. 4, Art. 7.13.2, formula 19). It may also be noted that, because of the assump- 
tion x > 2z, the first term in braces of (15) is negligible in comparison to the second. 
If in addition, x > n, then 


| T 2x \" 
(16) | Kn+ie(2) + K_n+iz(2) | sia x sinh (3x) (*) : 


From either (15) or (16), it can be seen that for m and z fixed, and x sufficiently 
large, an error of negligible magnitude can be obtained. For nm = 0, the error 
estimate is independent of z provided the condition x > 2 is satisfied. For ex- 
ample, with h = .5, we find that x = 4z, and an estimated error of 


1 
——— = 2.1 X 10. 
Ni sinh (47?) x 


| Kie(z)|~ 


In some problems it is more convenient to work with e*Ko(z) rather than with 
K,(z) itself. Evidently, to obtain this quantity to ‘Ye same accuracy will require 
a smaller interval if z is large. Thus h = .5 and z = .2 gives 


e*Ko(z) ~ 2.14075 73184 


as compared to the tabulated value from [8] of 2.14075 73163. If z = 10, an 
interval of A = .25 will give e*Ko(z) with a comparable error. By actual calculation 
the value .39163 1933 is obtained which agrees with Watson’s value [11] to as 
many places as are tabulated there. 

(4) The Error Function. This function is defined by 


His) = ~ f * Hat 


It can also be shown [6], that, 


y er ek 
H(s) = 1—-—e rs 





which is in a form suitable for evaluation by the present method. Taking 


—¢t? 


HO = sa 


then [4], 


g(x) = * e [=~ e-“dt f- a -at| ’ 
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For large values of x the first integral is nearly equal to f e-*dt = Vx while 
e7 (ettz)? 
22+ <x 


(2 + 4x*) 2 2x, 


. But since 





the second is of the order of 


therefore 
e~ > ete (etiz)? 


e7 (tHe)? 
22 + x 


edt 
P+? 





> 


with an interval “‘h’’ is of 





Thus the error in computing the values of f 
0 


the order of 


while the error for H(z) itself is approximately 


_2rs 


e *° 


Hence, to compute H(z) with a prescribed error of (say) e, the value of 4 must 
be so chosen that 


2x2 


e* <«e@ 


For example, if 22 = 2, we find that with h = .5, the error to be expected is 
e-"8 of 5 X 10-8, while with 4 = .25, the anticipated error will be 2.5 X 10-**. 
Actual calculations give: 


with h = .25, 
H(v2) = .0455002639. 


The second value agrees exactly with the one tabulated in [9]. 

In the present paper, we have considered only a few of the many functions for 
which the trapezoidal rule provides a simple and highly accurate means of 
obtaining numerical values. 

Henry E. Fettis 


Wright Air Development Center 
Wright-Patterson Air Force Base, Ohio 
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Trapezoidal Methods of Approximating Solutions 
of Differential Equations 


Introduction. Seemingly by historical accident the points of view usually 
adopted for stepwise methods of numerical solution of differential equations 
largely emphasize the discrete values found rather than the functions pieced 
together by the method over various intervals of advance. The authors, col- 
laborating as a professor and student team, have found that the straightforward 
selection of functions to approximate solution functions offers many advantages 
conceptually and makes it possible to see much-used methods in a new light. We 
present here a few of the simpler results. 

The purpose of this paper is to show that, by considering the method called 
the trapezoidal method (cf. MILNE [1 ] p. 24) as a parabolic or quadratic function 
method, not only does one obtain a satisfying geometrical picture of an approxi- 
mation curve in the usual case, but that two new trapezoidal or parabolic methods 
are suggested. Of the two new methods, one is based on a Gauss integration 
formula. Thus the new approach makes it possible to use the Gauss integration 
formulas. We believe that the two methods are new in their application although 
they are old in their use of polynomials. We have found in certain applications 
that both the new formulas have points of preference in some cases over the simple 
trapezoidal method. 

Quadratic Functions and Trapezoids. Let y’ = f(x, y) and let it be supposed 
that (xo, yo) is either an initial point or that we are now to treat it as such and 
that we wish to establish a value approximating the solution for x = x. Then, 
using Milne’s terminology, we write 


h 
(1) mu — y=s (y1’ + yo’), h = x1 — Xo. 


From (1) it is clear that, if the solution were known, then (1) is a trapezoidal rule 
for establishing the integral from the integrand values at the ends of the interval. 
In case of a differential equation in which f(x, y) involves y, however, we do not 
know 4; and if we agree to use y:’ = f(x, yi) we do not know y,’. Hence, (1) 
may be considered as an implicit equation for y; which, with adequate restrictions 
on h, yields directly an iteration procedure converging to a unique value 4. 
Now it is clear that in general 4;, as finally obtained, is mot exactly the value 
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of the desired solution function at x = x; and hence calling (1) a trapezoidal 
method is a concession to form rather than an exact interpretation of the trape- 
zoidal integration of a known function. To understand what is going on, we choose 
a function u(x) = do + ai(x — xo) + a2(x — xo)? which we want to use as an 
approximation to the solution from x» to x;. We require that 


u(xo) = Ye 
(2) u' (xo) = yo’ = f (xo, yo), 
u’ (x1) = f (x1, U1) 


where we write u; = u(x,). Geometrically, a quadratic function u(x) which 
satisfies these conditions gives a parabola through (xo, yo) which has the “‘direc- 
tion’’ f(xo, yo) at (xo, yo) and which has the field direction f(x:, u(x,)) at x = x. 
This function is interpreted as an approximation to the solution function for 
xo < x < x. The first two equations in (2) give ao = yo and a; = yo’. Using these 
in the last equation (2) we have an equation for a2 as follows: 


(3) 2ash = f(x1, Yo + yo'h + aeh*) — yo’. 


In general this equation cannot be explicitly solved for a2. However, if we as- 
sume that a solution a2 of (3) is determined and write f(x1, yo + yo'h + ah?) 
= f(x:, 41) = yi’ then we have on substitution from (3) for a2 


, , 

(4) u(x) = yo + yol(e — x0) +7" (@ — xa) 
With x = x;, u(x1) = yi we have exactly the “‘trapezoidal’”” method (1). We do 
not recommend using (3) to find a2; on the contrary, it is usually greatly preferable 
to use the iteration suggested by (1). However, in (4) we have, with the solution 
for y; and y;’, a quadratic function which may be considered as the approximating 
function from x» to x;. We have never seen this function mentioned in papers or 
treatises on numerical methods. Continuation of the solution method through 
several steps gives this picture: the approximating curve consists of continuously 
and smoothly joined parabolic arcs. Thus the method is theoretically symmetrical 
or reversible, a terminology we have never seen applied elsewhere. In some cases 
reverse integration of a differential equation offers some hope of estimating effects 
of roundoff errors, especially with theoretically reversible or symmetrical methods. 

While it is not necessary to interpret the trapezoidal method by a quadratic 
approximating function, it is simplest to do so. Now, as an integration formula 


(1) holds exactly for every quadratic function y(x). The error term has the form 
3 
— =" for continuously thrice differentiable functions y(x). We consider the 


quadratic approximation function 


(5) u(x) = yo + yo’ (x — Xo) + a2(x — Xo)’, 
and we ask in addition that 

(6) Up’ = f(Xp, Up), 

where 


Xp» =xo+ ph, up, = u(x,), and u,’ = u’(x,). 
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We consider the number # as a characteristic value to be chosen to our advantage 
later. Assuming a solution to (6) for u, = y, and letting u,’ = y,’, we then have 
that (5) and (6) give 


(7) u(x) = yo + yo! (x — x0) + are (x — x0)? 
Now letting u(x:) = y; we have 

h 
(8) = Not p Lye’ + (2p — 1)y0']. 


We require that (8) which automatically is satisfied for every quadratic shall hold 
for y = x*. But this then determines p = 2/3 or 


h 
(9) v1 = vot+ 4 (3y’2/3 + yo’). 


In another connection this may be recognized as HEUN’s formula used with 
y’-2, Yo’, y-2 and 4, respectively, replacing yo’, y’2;3, yo, and y:. Now, however, 
we make use of our quadratic function u(x). We cannot use (9) since y’2/; is 
not known but from (7) with x = xy, p = 2/3 we find, letting u(x») = yp, 


h 
(10) Yojs = Yo+ 3 (y’23 + yo’) 


which provides the same algorithm type as (1) for finding y’s;3; which is then sub- 


stituted in (9) to obtain ;. If (9) is considered as an integration formula it has 
hi IV 


an error term 316 assuming (x) has a continuous fourth derivative. The method 


is not symmetrical but it has been easy to calculate and more accurate than (1) 
in some calculations. 

The application of the Gauss formula results if we choose the quadratic 
function 


(11) u(x) = yo + ai(x — x0) + a2(x — xo)? 
subject to conditions: x, = xo + ph, xg = x0 + gh, 

Up’ = f(Xp, Up), 

Ue’ = f (xq, Ug), 


from which we obtain the function u(x), with u,’ = _ tbe’ = Yq’: 


(12) 


(13) u(x) = yo + De = £98 oe (x 


pb @ 
Again if we take x = x; in 98 we find with y,; = a 


x0) re - aA - Xo). 


(14) yn = Yo + —— p Lea - 1)yp’ + (1 — 2p)9,"]. 


~: 
Now we require (14) which holds for quadratic functions to hold‘also for y = x! 
and y = x‘. Asa result we obtain p = 1/2 — 1/12 and g = 1/2' + 1/12 (these 
are the zeros of the Legendre polynomial of degree 2 on [0, 1]}). Using the Jetters 
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rather than the numbers when convenient we find from (13) and (14) 


(15) Ye = Yot rrr [(2qg — P)¥p’ — Pye’) 

(16) Ya = Yo+ #.. [(2p — 9)¥e' — Gy’) 
2(p — q) 

(17) y= yo +5 Oe! +99). 


Hence, we use (15) and (16) to find by compound iteration y,’ and y,’ and then 
substitute in (17) to obtain y,. The iteration procedure given by (15) and (16) 


‘. rs] 
converges provided h < (v3 + 1)/m where m = max mal Note that, as inte- 


gration formulas, (9) and (17) are both trapezoidal but that the evaluation points 
do not coincide with the endpoints of the interval. The error in (17) as an integra- 
hey ® 

4320° 

The last method in continuation gives a picture of parabolic arcs joined con- 
tinuously but not smoothly, each of which agrees in direction with that indicated 
by the differential equation at two points which are symmetrically placed in the 
interval. Thus this method is symmetrical or reversible theoretically, while that 
of equations (9) and (10) is not symmetrical. 

Of the three trapezoidal methods which is best? The comparative study of 
methods for merit depends on so many features that the most one can hope for 
now is a qualitative judgment. The first method is simpler but not much more so 
than the second. If the symmetry is no consideration then one will expect the 
second method to be more accurate. The third method, on the other hand, is 
appreciably more difficult to use than the other two, although it is the most 
accurate and we believe it is also the most stable of the three methods. For 
automatic computing machines the use of irrational abscissas is no worse than 
use of, say, 1/3. We have tested these methods at the University of Wisconsin 
Numerical Analysis Laboratory and have found that for a linear differential 
equation the third method was significantly the most economical for a given 
accuracy. If roundoff error is likely to be bothersome, the presumably many 
fewer total number of steps required for the last method would minimize the 
difficulty. For general use on simple problems the second method should be better 
than the first. 

We have neglected to give methods for obtaining starting values for the 
iteration procedure. We have done so, since it seems that we can use (10) with 
y'2;2 = yo’ for the second method, and yo + phyo’ and yo + ghyo’ as first estimates 
of y, and y, for (15) and (16). The customary methods of guessing these values 
would make it more difficult to change h repeatedly, which these procedures may 
claim as an advantage. The high speed computing machines make it possible to 
enjoy the advantages of intricate methods. They may make it virtually necessary 
on difficult problems to use methods of increasing complexity in order to mini- 
mize errors. 





tion formula may be written 
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We will present later a geometrical view of certain classical methods, and other 
new methods, some of which are non-polynomial in character. This presentation 
is made to indicate the value of the general approach we have made. The parabolic 
methods are special cases of what we have termed the parametric methods which 
include certain methods of ApAMs, MILNE, MouLton, OBRECHKOFF, RiTz, and 
GALERKIN, for example. Extension of results of this paper to higher order inte- 
gration formulas is straightforward but would serve no useful purpose here as an 
illustration of our new point of view. 


The Advanced Numerical Analysis class of Mr. HAMMER has recently carried 
out rather extensive C.P.C. calculations to compare several methods for numerical 
solution of the differential equation y’ = x? + y? with initial point (0, 1). For this 
problem the third method here was superior to the simple trapezoidal method. It 
is intended to publish these calculations separately. Mr. ORVILLE MARLOWE 
carried out calculations on simple linear differential equations including y’ = y 
with initial point (0,1) and concluded that here the third method was most 
economical for a given accuracy, partly due to the fact that no iteration is needed 
for linear equations. 

PRESTON C. HAMMER 


Jack W. HoLLInGswortH 
University of Wisconsin 
Madison, Wisconsin 
General Electric Corporation 
Schenectady, New York 


Calculations in the University of Wisconsin Numerical Analysis Laboratory concerning 
methods i in this report were financed by funds of the Wisconsin Alumni Research Foundation. 
¥ 1W. E. MiLne, Numerical Solution of Differential Bquetions, John Wiley and Sons, New 
ork, 1953. 


Solving Systems of Linear Equations with a 
Positive Definite, Symmetric, but possibly 
Ill-conditioned Matrix 


Introduction. Often a system of linear equations to be solved has a matrix 
which is known in advance to be positive definite and symmetric. The normal 
equations for least squares fitting of a polynomial form such an example. However, 
if the polynomial is of reasonably high degree, the matrix of the normal equations 
is apt to be ill-conditioned. This may be seen by observing the origin of such 
matrices. In general they are of the following form: 


¥ xP ees HxPX7¥ xX 
i 


é 7 


Dx ++ OxZIT x2", - 
i i i 


Here the superscripts are exponents, N is the degree of the polynomial, and the 
x; are the values of the argument at which the data is given. For sums of high 
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powers the x; of largest magnitude tends to dominate and thus the last column 
tends toward a multiple of the next to last column. Thus the determinant ap- 
proaches zero as N increases, and the matrix becomes ill-conditioned. (See also 
HERZBERGER [1 ].) It is well known that for a system of equations with an ill- 
conditioned matrix, an erroneous solution can be obtained which seems to satisfy 
the system quite well. (See e.g., Saw [2], p. 23.) 

Various measures of the ill-conditioning of a matrix have been proposed. 
Perhaps the most common is the relative size of the determinant ; those with the 
smallest determinant being generally the most ill-conditioned. Other more precise 
measures have been proposed, which assume knowledge of either the eigenvalues 
or the inverse of the matrix. (See e.g., Taussky [3].) Among these are 





1 bee 
1 P(A) = 
©) 4) || min 
and 
(2) N(A) = ~ N(A)N(A~). 


Here A is the matrix of order n, \; its eigenvalues, and N(A) is the norm: 


(3) N(A) = i, a;;. 
The larger P(A) or N(A), the more ill-conditioned is A. 

The present paper describes a procedure for solving a system with a positive 
definite, symmetric, matrix, which (especially when used in conjunction with the 
Cholesky or square root method) involves little extra labor when the system is 
well-conditioned, forewarns the solver if the system is ill-conditioned, and extends 
somewhat the class of equations which can be satisfactorily solved without using 
double precision in the computations. A measure for the condition of the system 
of equations is proposed. Previous measures of condition have been only for the 
matrix itself. 

The Proposed Procedure. Let the system of equations to be solved be 


(4) AX = B, 


where A is assumed to be positive definite and symmetric. Let the eigenvalues 
of A be 
Amin =A SMS pint <= An = max 


with corresponding eigencolumns V;. The proposed procedure is to solve instead 
the perturbed system 


(5) (A+kI)Y=B 


with matrix C = A + kl, where k is a small positive constant, and then compute 
X from its series expansion in Y. The matrix C will be shown in Appendix I to 
be better conditioned than the matrix A. Present experience indicates that 10*-* 
and 10*-*, where s is the number of decimal places being carried, are reasonable 
values for R. 


fh Oe mes ta Se 
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Since A = C — kl, then 


(6) At=C1+ kCO?+--- +k ™C-™'14+.---; 
i.e., 


(7) X =A“B = CB+kC*B+--- + k"C™ B+... 
= V+kC1V4+--- + (RC)"V¥4+---. 


The eigenvalues of kKC— are k/ (A; + 2), and 0 < k/(A; + &) < 1 for all z and for 
k > 0, since all \; > 0. Thus the series converges. If Amin > R, then k/(A; + k) K 1 
for all 7, and the series converges rapidly. 

An Easily Obtained Indication of the Condition of the System. The procedure 
now is this: Using a value of k of one in the, say, next to last decimal place, Y is 
computed by any method, such as the Cholesky or square root method, which 
does most of the labor before the right member B is used. Then using Y as a new 
right member, kC—'Y is computed. This again may be used as the right member 
to compute the third term in the series. If at this point, it appears that further 
terms would contribute nothing in the range of decimal places being carried, then 
the implication is clear. The condition of the system must have been good, assum- 
ing k was appropriately chosen. In other words, the perturbation of the diagonal 
elements had small effect on the solution of the system. On the other hand, if 
each term remains relatively large, then the system must have been quite ill- 
conditioned, and probably no meaningful solution exists. If slow falling off of the 
terms is observed, then the system is probably fairly ill-conditioned. Direct 
application of Cholesky’s method to the system woyld probably produce an 
inaccurate solution. Even the accuracy of the Y computed is open to question. 


Use of the Method to Improve an Approximate Solution. It should be noted 
that, since 


) max + k 
P(A) 7 and P(C) = ha +h 


k must be approximately as large as \min or larger in order to appreciably improve 
the condition of C over that of A; while on the other hand, since k/(A; + R) are 
the eigenvalues which determine the convergence rate of the series, a k much 
larger than Amin will cause slow convergence. Thus only when A is fairly well- 
conditioned is k Y accurate enough to be used as a new right member to compute 
the second term in the series without computing instead B — AY, which equals 
kY by (5). 

However, there exists a class of matrices A which cannot be successfully 
inverted carrying a given number of decimal places, but for which the correspond- 
ing matrices C can be successfully inverted. The latter can be used to improve 
any approximate solution of (4). To see this, let Xo be an approximate solution 


of (4). Compute B — AX» accurately, and use this as a right member in the 
system. 


CZ) = B — AX». 
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Let X be the true solution of (4). Then 


CLX — (Xo + Zo) ] = CX — CXo — CZp = AX + RX — AXo — RXo — CZp 
= B+ kX — AXo — kXo — CZo = R(X — Xo) 
or 


xX —- (Xo + Zo) = kRC7(X > Xo). 


Since the eigenvalues of RC— are all less than one, if X; = Xo + Zo, and the 
process is continued, then X;— X. Since a significant inverse of C is available, 
an approximate solution can be improved without limit. (In this connection, see 
PoLACHEK [4 ].) 

If Xo = Y, the above processes are formally equivalent, but the second one 
does not have to contend with inaccurate right hand sides. 

In computing B — AX; more decimal places than previously carried are of 
course necessary. A desk machine is ideally suited to this purpose, as mulitplica- 
tion is not rounded off. However, those automatic calculators which can select 
the digits retained can also make this computation. Intermediate overflow will 
not in general invalidate the computation, since B — AX; is a residual vector 1 
and will be expected to be small. 

For accelerating the convergence of the sequence X; the 6?-process of AITKEN 
[5,6] and its extensions by SHAnks [7 ] are especially appropriate. (See For- 
SYTHE [8], p. 309-310). 

A Measure of the Condition of the System. If the elements of A are scaled 
for use by an automatic computer, then Amax has certain well-known bounds, so 
that P(A) can be large only if Amin is small. The significance of this is also seen if 
B is expanded in terms of the eigencolumns V;: 


(8) B= > bi Vi. 


i=1 j 


Then X = AB = D> 5b;V;/d;. Since B is rounded off to a fixed number, say s, 
i=1 

of decimal places, the b;V; are known to only s decimal places. If Amin = 10-*, 

then X can be determined to at most s — g decimal places. If the exact solution 

of the given approximate system has no significant digits in this range, then no 

significant solution is determined by the approximate system. 

Moreover, the vanishing of }, will not eliminate completely the effect of a 
small A; = Amin, Since small changes in the elements of A and B will cause changes 
in the Xx, Vi, and by. q 

Thus for a properly scaled set of coefficients, it is the size of Amin which is the 
essential criterion. 

The foregoing analysis suggests that P(A) is a satisfactory measure of the 
condition of the matrix A ; while on the other hand, the measure of the condition 
of the system of equations AX = B could be taken as, say 


P(A, B) = max > min =—— + P(A). 
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Here the eigencolumns V; are assumed to have been normalized. The purpose of 
the denominator is of course to remove the effect of scaling the system. 
Appendix I. Proof that C is better conditioned than A. 
1. det A < det C 


This is clear, since det A = [J \; and det C = | I (A; + &). 
i=1 #=1 
2. P(C) < P(A) 





) a k Amen 
This is also clear, since P(C) = wate and P(A) = we 
3. N(C) < N(A) 
The proof of this requires some of the properties of the norm given, e.g., in von 


NEUMANN and GoLpsTiNE [9]. For a positive definite, symmetric matrix 


N(A) = [x vf. 
Thus 


1 n ; aL n 1 ; 
N(C) = 1% 0+ 4) | [= aa | 


Na) ==| & wa} [Zs 


By Minkowski’s inequality (Harpy, LitrLEwoop, and P6LyYA [10], p. 31), 
for n > 1, 


and 


t=1 t#=1 


[= (as + we] < [= we + bn 








and 
n _ n -} 
[= (As + “| ‘a [= ne] + kn 
or 
° 1 4 1 
[= . | < n 5 
i=1 (Ay + k) | sl 4 kn- 
#=1 Ag 
Therefore 
: [= wef + kn 
N(C) < aT a 


n 1 -} a 
Loa] te> 
st 
n[ 5 Pl +k . n 
= es 7 < | x mp [ =| = N(A). 
m[ y "a + k t=1 s=1 AG 








i=t A? 


Appendix II. | Y| <|X|, where |X| is the length of the vector X. This is a 
property which is useful if the elements of the solution represent physical quan- 
tities which have definite fixed bounds. (See also LEVENBERG [11].) The proof 
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requires some properties of the upper bound of a matrix A, denoted by |A|, 
These properties are also available in von NEUMANN and GoLpsTINE [9]. 
For a definite matrix A 


|A| = Amaxe 
For any matrix A and a vector B 
|AB| <|A| |B}. 


Then, since 
Y= CB = (A + kI)B = (I + kA“)“A“B = (I + kA“)"X 
it follows that 


[¥|<|( + kA~)| |X| = max] —— 
5 dee 


ne ) are 
* max| <*é 5 |x! "Tg pect. 


James D. RILEY 





|X| 
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On the Numerical Solution of Elliptic 
Difference Equations 


1. Introduction. In recent years a considerable amount of progress has been 
made in improving rates of convergence of iterative methods for the solution of 
elliptic difference equations. The numerical solution of Laplace’s equation in 
rectangular coordinates in a unit square may be obtained by replacing the 
differential system by a difference system over a rectangular network with mesh 
spacing h, and then solving the resulting difference equations by an iterative 
method. This problem is usually taken as a “‘model problem”’ in elliptic difference 
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equations because the rate of convergence of iterative methods may often be 
obtained precisely. Prior to about 1950 the number of iterations for given accu- 
racy with various methods was about proportional to h~*. RicHARDSON [1] pro- 
posed a method with potentially greater rate of convergence, but gave no very 
effective way for choosing the parameters required. In 1950 YounG [2] obtained 
a method called the method of ‘‘successive overrelaxation’’ and independently 
FRANKEL [3] obtained his ‘‘extrapolated Liebmann’’ method, a special case of 
successive overrelaxation, and his ‘‘second-order Richardson” method, later 
shown by RiLey [4] to be also a special case of successive overrelaxation, if 
appropriate modifications were made. SHoRTLEY [5] obtained a method using 
Tchebycheff polynomials which was equivalent to Richardson’s method except 
that the parameters required were obtained by a definite procedure. Young [6] 
has analysed this method under more general assumptions about the matrix of 
the system. In all these methods the number of iterations for given accuracy is 
proportional to h-. 

It will be shown that further improvements in rate of convergence are pos- 
sible. Under certain conditions there are methods in which the number of itera- 
tions for given accuracy is proportional to A“ and even to h®. Specific examples 
of the use of the new methods will be given. 

2. General Description and Theory. Let our equations to be solved be 
C’V’ = p’. We normalize C’, transform our dependent variables, V’, and inhomo- 
geneous part, p’, by obtaining C = D“*C’D+, V = DtV’, and p = D-‘p’, where 
D is the matrix of diagonal elements of C’ and where the signs of our equations 
C’V’ = p’ have been chosen so that the elements of D are positive. We assume 
in all that follows that C is symmetric and positive definite, that C has what 
Young calls property (A) [7], and that the rows and columns of C have been 
ordered as in the extrapolated Liebmann method [3] (Young’s ordering o2 [8]). 
Self-adjoint linear second-order elliptic partial differential equations may always 
be replaced by difference equations in such a way that the preceding properties 
are obtained (Young [9 ]). 

The iterative sequence for the method of successive overrelaxation is 


t-1 N 
(1) Veta VM =a 2 eg¥if" +2 cgV™ — wl, 


j=1 i=i 
= 1,2,---,N;n>0). 


We define triangular matrices A and B such that A + B = aC and B = A’ 
— (2 — a)I. In terms of A and B, (1) becomes 


(2) AV+) + BV = ap. 


In (2) A is the matrix which multiplies the ‘“‘new” or “improved”’ variables in 
the extrapolated Liebmann method, B the “old” or “unimproved”’ variables. 
AT is the transpose of A. The difference, E™, between V and V™, i.e., the error 
in iteration n, satisfies 


(3) AE) + BE® = 0, 
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whence 
(4) E@+) = A-oBE®. 


If the length of E™ tends to zero as m tends to © we say that the transforma- 
tion A~'B is convergent [10]. The rate of convergence, Rx, is defined as — log K, 
where K is the spectral norm (magnitude of largest magnitude characteristic 
value) of A~'B [10]. Let i be the spectral norm of C — J, and let g = 1 — i. 
i lies in the range (0, 1) and may be evaluated for simple problems such as the 
“model problem’’ mentioned in the introduction. Young [2] has shown that: 
(a) K =|K,| where K, — a\Kii+a—1=0. (b) K =}? when a = 1. (We 
calla = 1 the Liebmann value of a.) If 7 is small Rx = 2. For the model problem 
with h small, ~@ = 342°h? and Rx = rh?. (c) K = K(min) = ao — 1, when 
a = ao = 2[1 — (1 — d*)#] + X*%. (We call ap the “optimum” value of a.) If i is 
small Rx + 2V2,, and, for the model problem, R = 2xh. (d) For a > ao, all the 
characteristic roots, K;, of A~'B have magnitude a — 1. 

The above results were also obtained by Frankel [3] for the model problem. 

So much for the ordinary extrapolated Liebmann method. We next consider 
the following iteration process: 


(5) AV@n+) + BV@™) = ap, 
(6) ATV nt?) 4+ BIV@nt) = ap, 


If we write (6) in ‘“‘component” form we obtain 


(7) V,&+2) = Vet) —a > Cij V;@**) op > Cij V,;@=+) — pil. 
j=i+l1 j=1 

In order to solve (7) for V;@"+®) we must first solve for Vy@"*”, then V9¥7?, and 
so on. Hence, if we think of (1) as representing the traversal of a network in a 
certain order (corresponding to the extrapolated Liebmann method) then (7) 
represents a traversal of the network in the : pposite order (the ordering, however, 
still corresponds to the extrapolated Liebmann method). Thus in the iteration 
process (5), (6) we are going back and forth over a network. We will call (5), 
(6) the method of symmetric successive overrelaxation. The error E®, in (5), 
(6) satisfies 


(8) AE@at+) + BE®») = 0, 
(9) ATE@™+) 4 BIEG@nt+) = 0, 
and hence 

(10) En) = AT BTA3ABE®), 


Let T = AT 'BTA-B. Then we have the following theorem : 
THEOREM 1. The characteristic roots of T are real, non-negative, and less than 
one for0 <a <2. 





104 NUMERICAL SOLUTION OF ELLIPTIC DIFFERENCE EQUATIONS 


Proof. 

T = AT’BTA“B 
= AT'[TA — (2 — a)IJA“[A? —- (2 — aI] 
= AT{[I — (2 — aA] - (2 - «ATA? 
= AT [I — (2 —a@)A“][ — (2 — a)A“)T}A?. 


The quantity in braces in (11) is a semi-definite matrix (because it is a matrix 
times its transpose) and T is the transform of this semi-definite matrix by a 
transformation that leaves its roots invariant. Hence, the roots of T are real and 
non-negative. We may easily show from Young’s results that K < 1 for0 < a < 2. 
(10) may be regarded as composed of two successive linear transformations each 
with spectral norm K. As two successive convergent transformations give a 
convergent transformation, T is convergent for 0 < a < 2, and hence the roots 
of T are less than one. 

The foregoing theorem is just a generalization to a in the range (0, 2) of a 
result obtained by AITKEN [11] for a = 1. The method of proof is similar to that 
used by Aitken. 

Let the characteristic roots of T be \;,7 = 1, 2, ---, W and let the spectral 
norm of T be \;. We have: 

THEOREM 2. If a 2 ao, \1 2 K?, equality obtaining if and only if the character- 
istic roots of T are all equal. 


(11) 


Proof. 
|A“B|-|A“B| = K24 
= |A-"BT|-|A—B| 
(12) = |T] . 
N 
= I ri. 


Therefore, Max \; 2 K?, equality occurring only if the \; are all equal, for other- 
wise we could not satisfy K?% = II,\;. The first equality in (12) follows from 
Young’s result (d) above. 

Coro.iary 1. If a 2 ao the rate of convergence of T is less than twice the rate 
of AB. 

Corollary 1 follows immediately from Theorem 2 and the definition of rate 
of convergence. Hence we have proved that per network traversal (5), (6) con- 
verges more slowly (except for the very special case where the }, are all equal) 
than (2), i.e., back and forth iteration is less efficient than going just one way. 

Let k = 1 — K and suppose that & is small. Then we expect that most of the 
roots \; will be grouped in a band about K°, for each root which is substantially 
less than K? requires many roots greater than K? in order that (12) may be 
satisfied. Hence we have the following rather loose corollary : 

CorROLLARY 2. If a > ao, if k is small, and if N is large, then almost all the ; 
lie in the neighborhood of K°. 

For the ordinary extrapolated Liebmann method, we find a simple relationship 
(Young’s result ‘‘a’”’ above) between K, the spectral norm of AB, and i, the 
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spectral norm of C — I. We would like to obtain a simple relationship between 
\, and \ for then we could make an exact comparison of the rate of convergence 
of the symmetric successive overrelaxation method with the ordinary successive 
overrelaxation method. However, using standard methods which work for the 
analysis of many other iteration methods, we find that we cannot even determine 
\, for our model problem (section 4). Hence it is doubtful that \, is related to \ 
through any simple equation. However, we can analyse the one-dimensional 
Laplace problem (section 3) and also we can perform experiments on a computer 
(section 3 and 4). This work leads us to the conclusion that the following assump- 
tion is plausible: 

AssuMPTION 1. If we let our mesh-spacing, h, tend to zero, k and wp; = 1 — ry 
tend to zero at the same rate, i.e., k and yu; are of the same order of magnitude. 

We will assume that our assumption is true. Then we can construct an intera- 
tion process which, for sufficiently small h, is superior to the extrapolated Lieb- 
mann method. This iteration process is 


(13) AVG) + BV) = ap 
(14) ATty™” + Bryent) = ap 
(15) Vent?) — yo + weal. yo —_ Ve). 


The w2, are numbers chosen according to a technique similar to that used by 
Shortley [5]. Let E& = V— Vem, FE = V— V"". Then E®, the error 
at iteration m, satisfies 


(16) AE@+) + BEG) = 0 
(17) ATE”™™” + BTEC) = 0 
(18) E2n+2) si E' Oe 4, w.(E™™* i E®), 


Expand E) in the characteristic vectors, U;, of the homogeneous linear system 
(16), (17): 


(19) Eon = ioe Ui. 
Then 
(20) C2) = [Ay + wam—n (Ai — 1) Jes" 
m—1 
= TIT Di + wn; — 1) Je. 
n=0 
Let 
m—1 
(21) Pa(d) = TL [A + wan(A — 1). 


n=0 


P,,(A) = 1 for } = 1. We may factor P,,(d) as follows: 


II @-A) 
(22) P(A) = =, 
II (1 — Aj) 


gut 
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where w2; = A;/(1 — A,;), as may be verified by substituting into (21). If nothing 
is known about the distribution of the \; except that they lie between 0 and \, 
then we let (cf. Shortley [5] and Young [6]) 


_ Tn(@r + 5) 
(23) Pp(d) = Tala +b)’ 
(24) Tn(x) = cos (mcosx), |x| <1 


cosh (m cosh—! x),|x| > 1. 


T(x) is Tchebycheff’s polynomial of order m. In the interval between ad + 0 


= — 1 and a\ + 5 = 1 the maximum value of P,,(A) will be 1/7,.(a + 5). We 








let ad + 6 = — 1 whendA = Oandad + b = 1 whend = Ay. We obtaina = 2/d,, 
21— 
b = — 1. The roots of T,,,(x) lie at x = cos y w,j = 1,2, ---,m. We obtain 
for A;, 
nN 2j3+1 
(25) As = © (cos Zt r+t). 
2 2m 


Let 1/T,,(a@ + 6) = r™. We call r the average spectral norm for the transformation 
involved in going from E®) to E@=+*), The average rate of convergence is 
R, = — log r. Suppose that 1 — r is small and that m is large. Then 





1 
(26) : o~oie [Tn (a ; ay | = — Flog n=) FI 
— hlog {e?yampium = Vu. 1 





We compute R,/2 instead of R, to get convergence rate per network traversal. 
Suppose that symmetric successive overrelaxation converged at the same rate 
per network traversal as ordinary successive overrelaxation. Then for our model 
problem we would have R,/2 = 2Vxh. We would not expect to obtain this rate 
but according to Assumption 1 we should have R,/2 proportional to hi. 

Young [12] shows that for the ordinary successive overrelaxation (with 
consistent ordering and a = ao) the error goes to zero as nK* not K*. The operator 
in (11) is diagonalizable so that in the new method the error goes to zero (on the 
average) as r™ not mr™ which gives the new method an additional advantage, 
especially when the accuracy to be obtained is not great. 

From equations (8), (9) we obtain: 


(27) AATE@*+?) — BTBEG) = 0. 
and for a characteristic vector, U; 


(28) (AAT; — BTB)U; = 0 


(28) is useful in analysing specific cases. 








(29) 


Equati 
(30) 
(31) 
(32) 
with | 
When 
(33) 
where 


(34) 


(35) 


(36) 


(37) 
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3. Laplace’s Equation in One Dimension. We have 


Vo = Ci 
(29) — $Vp1 t+ Vy — $Vo41 = 0 ¢- 

Vy = Ce 
Equations (13), (14), and (15) become: 
(30) vero om vy = a[} vet il vy + 1V,) = 0, 
(31) ae = ver ail af} vet if ver 4 Ve = 0, 
(32) ver? = V5" + wl V5” — V5"), 


with Vo and Vy equal to c; and cz wherever they occur in (30), (31), and (32). 
When we work out equations (28) for this problem we obtain: 


(33) O5-1U 5, 5-1 + b,U i.» + 0nU i; p41 - 0, 
where 


(34) ao = ay-1 = 0, 


(35) a a,= —5u+5@-1), = 1,2,---,N —2, 


(36) b=), 


a a 
(1+£)ax-(14+2) +00 - a, p =2,3,---,-N—1, 


(37) b =A; —1+ a(2 —a). 


We try as a solution to (33): 


(38) Us» = sin— 





Bert. 


All the equations (33) have the same coefficients except the first and the last. 
If any one of the intermediary equations is satisfied by (38) then all will be and 
so will the last equation. We will also have: 


(39) aU.o + bUi1 + aUi2 = 0 
which may be combined with the first of equations (33) to obtain: 
(40) aUu;.o + (b _ b:) U1 = 0. 
We substitute (38) into the intermediary equations (33) and obtain for \,, 
1+ 2 —a(2 — a) — (a — I)at; 


(41) i - 2 ’ 
1+7 — ats 





in : , . , ade 
where ¢t; = cos W B;. It is convenient to introduce the following quantities: 
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(42) wm=1-X; 
(43) w=1-44, 
(44) 6=2-a 
82 
(45) &; = rowel 


With (42), (43), (44), and (45) in (41) we obtain 








6 
46 i= . 
(46) MITE 
The boundary condition (40) becomes 
N-1 
(47) t; sin Bri + : sin —— Beri = 0. 


We choose 6 so that &; = 1. We can show a posteriori that this gives closely 
the optimum value of a when N is large for both (5), (6), and (2). We shall 
assume WN large and obtain an asymptotic solution for the y;. First we consider 
when 7 is small. Let 8; = 1 + «;. Then from (47) 


a a 
4 aa se T =, * 0. 
(48) (S+a)ats . 
Weh pe ae 
e a a so ' 
6? 7¥1 1 
49 {= — & -, 
(49) : 4ayi Yi e 
1 
(50) o  —— 


~ a/2 + & 2N 
Hence ¢; is small for the first few values of i where the expansions used are valid. 
6 , — 
We have pi = 3 and the successive values of wu; decrease as 7 increases, and 


é; rapidly grows small. Next we obtain a solution when &£; is small. Let 


N 
8; = Wn (1 + «;). Then 





(5 1) &; sin 





a Fr wi(1 +e) + = sin ri(1 + €;) = 0 


(52) sin wi(1 + ¢;) = (— 1)*wie; 





N 
(— 1)*€; sin Vu mt 


(53) ¢ = ; , +#N—1. 
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N 
W for the larger values of i and so y; — 2, &; — 0, 
and the y; approximate 6 for large 7 so that we have a narrow band of roots from 


6 
—gtoi-—-- 
1 oO 2 


Hence 8; is close to 


The approximation in (51) of neglecting ¢; in the first term to obtain (53) is 
not valid when i = N — 1, for then we have 


(54) tv_1 sin [Nw (1 + ev_-1)] + : sin [(N — 1)4(1 + ev-1)] = 0 


and there is no real solution with |ey_,| < 1, when &y_, is small. The last charac- 
teristic value is obtained from 


(55) Un-1,p = sinh (N - P)Bn-1. 
For this characteristic vector, ty_1 = cosh By_; and yy_1 is negative. (40) becomes 
(56) 6? sinh NBy-1 + 2a?(1 — cosh By_;) sinh (N — 1)By_1 = 0 


and if By_; is small, 


(57) (8? — a’6*y_,) sinh NBy-1 + a’*y_; cosh NBy_; = 0. 
Let 

82 
(58) Pun = — (1 + €). 

a 


Then from (57), (58) 


6 cosh x 6 
(59) «= = — 


asinhr a 





and when N is large so that 6 is small we can show that 
(60) gy... = 1 — 6 
and so uy-1 = 1 and there is a single root close to zero. 


R, R, 
We obtain for ry R,, and * + R,: 





R, 
(61) tn Vah 
(62) R, = 2xh 
R, 1 
= -—— h-. 
(63) 2R, 2vx 


The factor of improvement is not large unless h is quite small, N quite large. 
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For this simple one-dimensional problem we can construct an iteration process 
in which the rate of convergence is independent of N. If we choose, for example, 


a 
w, alternately equal to zero and — 1, then the first w will damp out the error 
with characteristic value near zero and the second will damp every term in the 
6 , 
band (1 -—6;1- *) by at least 3. Further by choosing a smaller value of £, at 


the start, i.e., a larger value of a, we may increase the rate of convergence without 
limit, but the w, also increase without limit so that the whole process becomes 
unstable because of rounding error. 

An example of this iteration process was tried on the IBM 701 computer [13]. 
N was chosen to be 100 and c; and ¢2 were chosen to be zero. A starting approxi- 
mation of V equal to 10° at all interior points was chosen. The w, were chosen 
according to the scheme suggested in the preceding paragraph. Using the ordinary 
extrapolated Liebmann method approximately 300 iterations were required ito 
damp the maximum error by 10-5. Using the new method 16 iterations (32 passes 
over the network) damped the maximum error by about the same amount. 

By using a polynomial derived from the Tchebycheff polynomials, as in sec- 
tion 2, to damp the amplitudes of the error characteristic vectors with roots in the 


) : ’ 
band (1 —6;1-— ‘) we could obtain an even greater factor of improvement. 


We make this analysis of Laplace’s equation in one dimension in order to 
learn about the properties of symmetric successive overrelaxation, but not with 
the expectation that this iteration method would actually be used for one-dimen- 
sional problems. For one-dimensional problems direct Gaussian elimination 
methods are most efficient [14]. ‘ 

4. Laplace’s Equation in Two Dimensions. A detailed analysis of the model 
problem mentioned in the introduction does not appear possible. We may obtain 
equations (28) in explicit form for the two-dimensional problem without much 
difficulty. We find that we can construct vectors composed of trigonometric and 
exponential functions which satisfy (28) at all mesh points not adjacent to 
boundaries, but then we find that these trigonometric and exponential functions 
may not be chosen in such a way that the \, are real and the boundary conditions 
are satisfied. 

The two-dimensional problem was investigated by solving a numerical ex- 
ample on a computer [13]. A 30 X 30 network (including the boundary points) 
was chosen for the problem of the square for which the optimum value of a is 
a = 1.805. A starting approximation of V“ = 10° was chosen at all interior 
points and the boundary condition V = 0 was applied over the entire boundary. 

The ordinary extrapolated Liebmann method required 110 iterations at 
optimum a to damp the maximum error by 2.5 X 10~-*. By iterating with the 
process (1), (2) (i.e., w, = 0) it was determined that \, = .88. A set of w, for 
m = 20 was constructed by the method suggested in section 2. Then the process 
(1), (2), (3) was applied for 20 iterations (40 passes over the network). At the 
end, the maximum error was damped by approximately 2.5 X 10-®. This reduc- 
tion in error corresponded quite closely with the theoretical value. Hence we 
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find that we have an improvement of a factor of about three for this two-dimen- 
sional example. A factor of three can be very important when a long series of 
problems is to be solved. However, we note that the method of symmetric succes- 
sive overrelaxation requires storage for two iterates, V’"” and V®”), whereas 
ordinary successive overrelaxation only requires storage for one. 

A distribution of w, similar to that chosen for the one-dimensional problem 
was tried for the two-dimensional problem. Rapid convergence was not obtained 
indicating that there are roots more or less evenly distributed between \, and 0. 
However the results obtained lead us to conjecture that there is an iteration proc- 
ess of order h® even for the two dimensional problem. Label the corners of our 
square a, b, c, d in a clockwise direction. Suppose that our equations are ordered 
so that the first equation is centered at the meshpoint just interior to a, the last 
just interior to c. When we apply a distribution of the o, like that of section 3, 
we find that we have rapid convergence at mesh-points along the line ac but slow 
convergence (or even divergence) in the neighborhood of the corners 6 and d. 
Let a — c indicate that we pass over the network using the extrapolated Liebmann 
method starting at a and ending at c, and similarly for other letter combinations. 
(There are two ways of doing this but it should not make much difference which 
we choose so long as we are consistent.) Starting with V™ we obtain V“™ by 
iterations a > c, c >a, b +d, d + b. Then we obtain V“™ from 


yor) — yer wr w,(V""" — vm), 


It seems plausible that a choice of the w, analogous to what was done for the 
one-dimensional problem would give convergence of order zero for this process. 

5. Suggestions for the Use of the New Method for Multi-dimensional Prob- 
lems. It is supposed that the method here presented would be used mostly in 
situations where we wish to solve a number of similar problems. In this case the 
following procedure is suggested : 


(1) For a typical problem obtain a and \, by experimentation on a computer. 

(2) Estimate \,(max) for the series of problems. 

(3) Determine m and the w; from \,(max) so that r™ is less than the desired 
accuracy e to be obtained. 

(4) Use the values of the w; obtained for the remaining problems. 


\, (max) should be chosen conservatively so that it is sure to be large enough 
to cover all the cases treated. If \; is chosen too small the new method becomes 
inefficient, so that it is better to overestimate \,. It is not believed that the new 
method requires a very precise determination of a. 

The new method is especially suited to a computer which has magnetic tapes 
which read “‘backward”’ as well as “forward.” 

The new method will perhaps also be useful for the numerical solution of 
elliptic difference equations in cases where relaxation methods converge very 
slowly ; for example, for the equations of compressible fluid flow at Mach numbers 
close to but less than one. 

Joun W. SHELDON 
Wilson Point 
South Norwalk. Conn. 
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Solutions of a Nonlinear Differential Equation 
Arising in the Theory of Diffusion Flames 


An extended theory of diffusion flames [1] leads to a nonlinear second order 
differential equation, which is, in its non-dimensional form [2 ]: 





v”’ : -v-(v+x+ x, — 1) 


etme 
with the following two-point boundary conditions: 


v=i1-—x, at *«=0 
and 
v—-0 as x—-o, 


This equation describes the reaction of sodium vapor diffusing into a halide, where 
1 — x; and K correspond effectively to a measure of the nozzle opening and the 
chemical rate of reaction, respectively. The equation has been solved numerically, 
with the aid of FERUT, to cover a wide range of experimental conditions. It 
presented many difficulties of a practical nature, which could hardly have been 
foreseen until considerable work had been done. 

First of all, solutions for two cases, corresponding directly to specific experi- 
ments, were carried out by hand, using the equation in its original form [1] and 
employing M1Lne’s method of integration. Each of these required nearly a month 
of desk work by a proficient computer. One of them served as an independent 
check for the FERUT programs. 

In order to adapt the problem for automatic solution it was decided to work 
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over a mesh of parametric values of K and x;. Correspondence with experiment 
could then be obtained by interpolation into these results. 

In an attempt to draw up a comprehensive program, certain difficulties soon 
became evident. For instance, because of the large variation in the parameters, 
the functions involved extended over a wide range of magnitudes, and in a manner 


x | 


Approximate relationship : 
1. Log Ka logS for fixed x,. 
2. Separation between curves o& log (1-2c)) 
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which could only be predicted from numerous and prohibitively tedious hand 
calculations. Also the effective value of « was not known, but soon appeared to 
be significantly dependent on the combination of parameters used. This rendered 
impracticable the technique of integrating both outward from x = 0 and inward 
from x large. It was also noted that the solution was only defined for a particular 
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value of the initial slope v’ (0). This necessitated an iterative or repetitive numeri- 
cal method, which at once introduced questions of convergence. 

It was decided to use the Runge-Kutta process of integration as modified by 
Git [3], a routine for this already being available in library form. This method 
does not require the separate calculation of starting values and cuts inherited 
errors to a minimum. However, it is slow in comparison with other methods [5]. 
The technique of repeated outwards integration, with systematic adjustment of 
initial slope to satisfy conditions at ©, was adopted. 

A fixed-point program, supplemented by a plotting routine to display on 
paper the shape of the approximating solutions, was first prepared. The method 
of adjustable pre-scaling [4] was employed. Scale factors were provided for x, 
and K, and it was hoped that these could always be selected to prevent numerical 
overflow and still give the desired accuracy. This program was used with some 
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success for the smaller values of 1 — x; and K, but proved inadequate elsewhere. 
Either overflow occurred regardless of choice of scaling factors or the interval of 
integration could not be suitably refined. Sufficient work was done to reveal that 
the entire solution, and not just its terminal behavior, was very sensitive to slight 
changes in the initial slope. This effect became progressively more pronounced 
towards the larger values of K and 1 — x; and made it necessary to recommence 
each integration at x = 0. Also, the converged value of the initial slope varied 
considerably from one parametric case to the next. Numerical details to supple- 
ment these various statements are summarized below. 

In view of these many practical difficulties, it was decided to make a direct 
attack on the problem with a floating-point program. This was actually written 
in Transcode, a system of automatic coding recently developed for FERUT. It 
paralleled the fixed-point program with the exclusion of the plotting routine, and 
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proved much easier to develop and operate. The reduction in running speed by a 
factor six was not intolerable. 

Some operating details follow : The initial slope, v’ (0), was defined as —s, since 
it was always found to be negative. Preliminary runs for each parametric case 
were made as follows: 


1. To determine the interval size which would yield the value of v(x) to 
between 3 and 5 significant figures. 

2. To determine two values of the trial s, one producing divergence to + © and 
the other divergence to — ~. 

3. To determine the effective value of x — © for applying the test for di- 
vergence. 


Final runs for each parametric case were made to determine s such that (x) 
dropped to within 1/200 of its initial value and remained parallel to the x axis 
for at least 1,20 of the total range of integration. The value of x at which v(x) 
first reached 1/200 of its initial value was defined as the effective value of «. 
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In the fixed-point program, s was automatically adjusted by averaging the two 
most recent trial values known to produce divergence in either direction. This 
proved rather crude in view of the sensitivity of the solution to the trial value of s. 
With the floating-point program, adjustment of s was made by higher order 
interpolation away from the computer. Since each run over the range of integra- 
tion required from 5 to 15 minutes, this procedure could be followed quite effi- 
ciently. Preliminary runs required integration over about 20 to 50 steps. However, 
as many as 180 steps were necessary in some of the final runs, when the value of s 
was approaching convergence. Also, the interval size could not be varied suitably 
over the range in view of the extremely divergent behavior of the approximating 
solutions. 

Figures 1 and 2 illustrate the relationships discovered between the converged 
sand each of K and 1 — x. Figure 3 presents the graphs obtained for 1 — x; = 0.2 
and the various K. Graphs for other values of 1 — x, display the same general 
forms. Graphs and tables of v(x) are available for all the following cases: 


— x, = 0.01 with K = 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, and 5. 
— x, = 0.2 with K = 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, and 5. 
= 0.4 with K = 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, and 10. 

— x, = 0.6 with K = 0.01, 0.05, 0.1, 0.5, 1, 5, and 10. 

— x, =0.8 with K = 0.05, 0.1, 0.5, 1, 5, 10, 50, and 100. 

— x, = 0.9 with K = 0.1, 0.5, 1, 5, 10, 50, and 100. 


ee 
| 
& 
= 


The experimental test case corresponded to 1 — x; = 0.004970 and K = 0.03713. 

Initial slopes ranged from —0.000296 (for 1 — x; = 0.01 and K = 0.001) to 
—3.463 (for 1 — x; = 0.8 and K = 100). Intervals of integration necessary to 
give 3 decimal accuracy ranged from 1.0 for the smaller 1 — x; and K, to 0.008 
for the larger values. The effective value of x — , as defined above, ranged 
similarly from 173 to 1.25. A reversal of the general trends was noted as s became 
greater than unity. 

The author is greatly indebted to Professor C. C. GorTiieB for advice in 
preparing this manuscript, and to Mrs. MAry BurGEss, who organized and carried 
through the production schedule for the floating-point program, and displayed 
much ingenuity in performing the interpolations. The chemical significance of 
the solutions is to be discussed in a separate paper [2]. The problem was sub- 
mitted to the Computation Centre by the National Research Council of Canada. 

BEATRICE H. WoRSLEY 
Computation Centre 


University of Toronto 
Canada 


1 R, J. Cveranovic & D. J. Leroy, ‘“‘Contribution to the theory of diffusion flames,” Canadian 
Jn. Chem., v. 29, 1951, p. 597-603. 


2R,. J. Cvetanovié, to be published. 
_ ,*S. GILL, “A process for the step-by-step integration of differential equations in an automatic 
digital computing machine,’’ Cambridge Phil. Soc., Proc., v. 47, 1950, p. 96-108. 
4B. H. Worstey, ‘Numerical representation in fixed-point computers,’ Computers and 
Automation, v. 4, No. 5, 1955, p. 10-13. 
5 In retrospect, it would appear that a MILNE-type formula, with programmed calculation of 
starting values, could have been more profitably employed. 
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Remark on Determination of Characteristic Roots by Iteration 


In this journal [1] one example was given of a 4th order symmetric matrix 
for which a certain algorithm did not quickly or conveniently furnish the values 
of the characteristic roots; indeed one of the trial values converged to a negative 
limit by “‘passing through + 0.” 

The same matrix was subjected to an established machine routine in this 
laboratory for the determination of characteristic roots and corresponding charac- 
teristic vectors by iterated orthogonal transformation. The computing procedure 
is due essentially to Jacobi [3] and is very similar to that reported in this journal 
[2] by Rospert T. Grecory of the University of Illinois; it was called to our 
attention by Dr. H. H. Go.pstine of the Institute for Advanced Study. 

The computational procedure employed is that of applying sets of successive 
plane rotations, each of which reduces to zero one of the non-diagonal elements 
of the matrix, until the sum of squares of non-diagonal elements fails to decrease. 
For a matrix of order n, the number of plane rotations in each set is $n(n — 1), 
one for each of the non-diagonal element positions. The computation for each of 
the rotations requires multiplication by two matrices which differ in only four 
places from the identity matrix, and, because of the symmetric nature of the 
original matrix, the number of multiplications needed for each rotation is only 4n. 

The matrix in question was 


2 1 3 4 
i =-3 1 5 
a= 3 1 6-2 
4 > =-2 -!1 


and five sets of six rotations each, i.e., not more than 30 plane rotations, pro- 
duced the eigenvalues 


Ai = +5.6688643, A, = —1.5731907, As = +7.9329047, dy = —8.0285783; 


an independent computation showed that the error in each value is less than 
1 X 10-7. A simultaneously performed auxiliary computation involving 4n addi- 
tional multiplications gave the eigenvectors as columns of the matrix 


+.3787 0268 —.6880 4793 +.5601 4450 —.2634 6239 
+.3624 1904 +.6241 2285 +.2116 3276 —.6590 4071 
—.5379 3516 +.2598 0086 +.7767 0826 +.1996 3352 
+.6601 9881 +.2637 5026 +.1953 8161 +.6755 7335 


An independent (hand) computation verified that the matrix product 
U (Diag. (+.5.6688643, —1.5731907, +-7.9329047, —8.0285783))U? = A, 


differed from A in all sixteen elements by less than 4-10-’. 
Similar (hand) computation evaluated the characteristic polynomial f(A) 
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= \* — 43 — 73d? + 260A + 568 in the neighborhood of the stated roots and 
showed that the error in the eigenvalues was less than 1 x 10-’. 
The last component of x in [1] should read —0.38333123+ in place of 
—0.3833123+-. 
J. L. BRENNER 


G. W. REITWIESNER 
Computing Laboratory 
Ballistic Research Laboratories 
Ordnance Corps, U. S. Army 
Aberdeen Proving Ground, Maryland 


Sponsored by the Office of Ordnance Research, U. S. Army 


1E. BopEwie, “A practical refutation of the iteration method for the algebraic eigenproblem,” 
MTAC, v. 8, 1954, p. 237-239. 
*R. 


. GREGORY, ‘ ‘Computing + wor and eigenvectors of a symmetric matrix on the 
wees 2 ¥ “MTAC, v. 7, 1953, p. 215-22 


ec 2 Jacosr, “Ein leichtes not .. .”, Jn. reine angew. Math., v. 30, 1846, p. 51-95. 


A Note on the Summation of Chebyshev Series 


In a recent note [1 ] I gave tables of the numerical coefficients in the Chebyshev 
expansions of some common functions, and suggested that these should be used 
to provide polynomial approximations by truncation and rearrangement in 
powers of the independent variable. 

The purpose of the present note is to show that the truncated Chebyshev 
series may be evaluated directly without rearrangement, and without using tables 
of Chebyshev polynomials. The process used is one of recurrence, and is well 
suited for use with automatic computing machines. 

Let the truncated series be denoted by 


(1) f(x) = Ao + AiTy* (x) + AsT2*(x) + --- + AnTy*(x), 
where, following the notation of Lanczos [2] which is used throughout, 
T,*(x) = cosn0, cos@=2x—1, O<x<1l. 


Ss 


In order to evaluate f(x) for a given value of x we construct the sequence By, 
By-..1, «++, Bo, given by Bys2 = Byy; = 0 and 


(2) B, — (4% — 2)Basi + Bayz = An, 1» = N,N —1,---, 1,0. 
Then we have 
(3) f(x) = Bo — (2x — 1)B,. 
This result may be verified by using the recurrence relation 
T.*(%) — (4 — 2)T*nyi(x) + T*ny2(x) = 0. 


It may be thought that the coefficient (4x — 2) in (2) can give rise to con- 


siderable building-up error in f(x). To investigate this point, we observe that the 
general solution of the equation 


Un — (4x 2 2) tensa + Unze = 0 
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is given by 


u, = aT,*(x) + BU,*(x), 


where U,*(x) = sin (m + 1)@cosec 0, cos@ = 2x — 1, and a, 8 are arbitrary 
constants. 

Hence a rounding error e() in A, or B, introduces an error e,(m) in B,(s < n), 
given by 


(4) é.(n) = aT ,*(x) + BU,*(x), 


where 
e(n) = aT,,*(x) + BU,*(x) \. 
0 = aT* n41(x) + BU* ns1(x) 


Solving the last pair of equations, we obtain 
a = e(m)-U*n41(x)/Ti*(x), B = — e(m)-T*n41(x)/T1* (x) 
and substituting in (4) we find that 


(5) és(m) = {T.*(x)-U*ngi(x) — Us* (x) T* ngs (x) }e(m)/Ti* (x) 
= U*,_,(x)-e(n). 


The corresponding error in f(x) may be found from (3), and is given by 


eo(m) — (2x — 1)ex(m) = {Un*(x) — Tr*(x)- U*s1(x) }e(n) 
= T,*(x)-e(n), 


which is the same as that produced by an error of e(m) in A, when the series is 
summed in the usual way. Hence the rounding-off of B, to the same number of 
decimals as A, can only double the maximum rounding error in f(x). If one or two 
guarding figures are retained in A, and B,, this error may be made negligible 
compared with the truncation error which has been introduced in replacing the 
infinite Chebyshev series by f(x). 

Thus, although the errors in B, may become quite large if N is large (see (5)), 
the error in f(x) is not serious. 

We may note that a similar device may be used to evaluate other series in 
terms of functions which satisfy a linear recurrence relation; for example, Neu- 
mann series of Bessel functions. 

Let 

f(x) = aopo(x) + aipi(x) + --- + anwpn(x), 


where 
(6) Pa+i(x) + anPn (x) + BaPn—i(x) = 0, 


and ap, 8, may be functions of x as well as of n. 
We construct the sequence by, by-1, ---, bo, where by41 = byy2 = 0 and 


bn + aadasi + Bazibdnye = Gn, n=WN,N— 1,---, 0. 
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Then application of the recurrence relation (6) shows that 
F(x) = bopo(x) + bif{p1(x) + aopo(x)}. 


C. W. CLENSHAW 
National Physical Laboratory 
Teddington, Middlesex 
England 


This note is published with the permission of the Director of the National Physical Laboratory. 
1C, 7: CLENSHAW, “Polynomial approximations to elementary functions,’ MTAC, v. 8, 1954, 
p. 143-147. 


2 NBS Applied Mathematics Series 9, Tables of Chebyshev Polynomials S,(x) and C,(x). U. S. 
Govt. Printing Office, Washington, 1952. 


Conjectures Concerning the Mersenne Numbers 


Conjectures concerning the Mersenne numbers are appropriate since they 
were inaugurated with one. A conjecture [1] that seems likely to be false, but 
unlikely to be proved false, is that all numbers p, are prime (n = 1, 2, 3, ---), 
where, for example, p, is 


Recursively, ~; = 3, Pay; = 29" — 1. The first four are 3, 7, 127 and 27 — 1, 
all known to be prime. Any factor of ps is congruent to 1 modulo 4, so ps cer- 
tainly has no factor less than 2!27. Similarly 


1 — 1 


is not divisible by any known prime, if 2?°*! — 1 is still the largest known prime 
[2]. One can try to argue about the probability that a number of the form 
2? — 1 is prime, when # is known to be prime. The probability that a whole 
number x is prime is about 1/log x, and is close to 


(1) jd — Ha - Na -H--- (1-4) 
where g = yx, so the factors (1 — 4), (1 — 4), etc., can be regarded as prob- 
abilities that are not far from independent. But if x = 2? — 1, only every pth 
factor of (1) should be taken, and the probability apparently ought to be about 
the pth root of 1/p-log 2, which is approximately 1 when ? is large. But this 
argument is also invalid, as we may see from the statistics of Mersenne primes 
[2]. We may see from these statistics (assuming them to contain no gaps), that, 
if m, denotes the mth Mersenne prime (m, = 3), then 


2.18 log log m, < m < 2.72loglogm, (3 < m < 17) 


while 


2.31 log log my; = 17. 
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It is reasonable to suppose that the number of Mersenne primes less than <x, 
when x is large, is about 2.3 log log x. This conjecture may be shown to be equiva- 
lent to the assertion that the probability of 2? — 1 being prime, when # is known 
to be prime and is large, is about 1.6(log p)/p, and is perhaps asymptotically 
(loge p)/p. If so, the probability that ; is prime is negligible, and we should be 
able to say with confidence that our original conjecture was the exact opposite 
of the truth. 

I. J. Goop 
25 Scott House 
Princess Elizabeth Way 
Cheltenham, England 


1E. CaTaLan, Nonu. Corresp. Math., v. 2, 1876, p. 96; cf. L. E. Dickson, History of the theory 
of numbers, v. 1, 1934, p. 22, ref. 116. 
?:—D. H. Lenmer, MTAC, v. 7, 1953, p. 72. 
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55(A, F].—Horace S. Unter, “Hamartiexéresis as applied to tables involving 
logarithms,”’ Nat. Acad. Sci., Proc., v. 40, 1954, p. 728-731 [1]. 


Hamartiexéresis appears to be a technical term in theology, meaning the 
absolute removal of sin. 

This paper contains in tabular form, the exponents of the prime factors 
(2, 3, ---, 997) in the product (1!) (2!)---(1000!). 

This table was used to check the first thousand entries in the table of F. J. 
Duarte [2 ]. Two errors were found: 

log 99!: the seventh quartet should read 8029 instead of 8929. 

log 266!: the eighth quartet should read 1897 instead of 1987. 

Later calculations indicate no (non-cancelling) errors in the range from 
n = 1001 to m = 1200. 


a2. 
1 See also Nat. Acad. Sci., Proc., v. 41, 1955, p. 183, for errata. 
2F. J. Duarte, Nouvelles tables de log m! 4 33 décimales, depuis m = 1 jusqu’A m = 3000. 
Geneva and Paris, 1927. 


56[C, D, E, K, L, S ].—Crcit HastinGs, JR., JEANNE T. Haywarp, & JAmes P. 
Won, JR. Approximations for Digital Computers. Princeton University Press, 
Princeton, N. J., 1955, viii + 201 p., 25 cm. Price $4.00. 


This book contains rational approximations of the following functions with 
approximate precision as indicated (there are several approximations to each 
function and the approximate precision of each is shown) : 


logio x, 10+ < x < 10, 3D, 5D, 6D, 7D; g(x) = (1 — e*)/x, 0S x < @? 
3D, 4D, 5D; arctanx, —’'1 <x <1, 3D, 4D, 5D, 6D, 7D, 8D; sin 42x" 
—1<x<1,4S,6S,8S;107,0 < x < 1, 4S, 6S, 7S,9S; W(x) = e*/(1 + e-*)* 
— «0 <x < ©,3D,4D,5D; E(x) = e*/V2r, — © <x < »,3D,3D,4D; 
K(n) = (n — 2n® — 2n*) In (1 + 2/m) + (2m + 18n? + 16n* + 4n*)(2 + m)~, 
O<n<o, 3D; T(i+x), O0< x <1, 5D, 5D, 6D, 7D; ¥(x) = (4/2 
— arcsin x)(1 — x)+,0 < x < 1, 4D, 5D, 6D, 7D, 8D; log. x, 2-4 < x < 23, 
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2 Pr 
5D, 8D, 10D; @(x) = =f e-*dt,0 < x < »,5D, 6D, 7D, and by another 
T 


formula 4D, 4D, 7D; K(k »¢ Sk <1, 5D, 6D, 8D; 


fF: dey 
do Vi — # sin? ¢ 


a/2 
E(k) -f Vi — R sin? ody, 0 < k <1,5D,6D,8D;In (1 +x),0<x <1, 


2 ot 


4D, 5D, 6D, 7D, 8D;e*,0 <x < ~,4D,5D,6D,7D; —Ei(—x) -f > 


z 





1 ” ‘. 
1<x < o, 4S, 6S, 8S; x(q) where g = mM edt, 0<2<.5, 3D, 4D; 
T z(q@ 
sad * du “sin (t — x) 
_ —,0<2<0,4D,5D;P -{ renee 
ws) i} ki(u)+rie(u)u’ ~ piner ewes 


1 <x < ,4D,5D,7D;Q(x) = f oe < x < ©,4D,5D, 7D. 

The notation above is that of the authors. In particular their E’ (x) is denoted 
as 2(x) in FLETCHER, MILLER, and ROSENHEAD, An Index to Mathematical Tables. 
The authors’ ®(x) is H(x) in the Index. The notation with regard to precision is 
an adaptation (by the reviewer) of notation used for precision of tables ; kD means 
that the absolute value of the difference between function and approximation 
does not exceed 5 X 10-*, and RS means that the absolute value of the ratio of this 
difference to the value of the function does not exceed 5 X 10-*. 

The function K(m) is involved in KLEtnN-NISHINA cross section determinations 
for nuclear scattering. The function W(z) occurs in some aerodynamical prob- 
lems [1 ]. 

Many of these approximations have appeared previously in MTAC, v. 7 
and v. 8. 

The approximations are in a sense equivalent to tables with a minimal number 
of entries and complicated rational interpolation formulas. The tabular entries 
are so few in number that all must be used in any interpolation ; these approxima- 
tions are equivalent to the resulting interpolation formula. 

The formulas are developed by fitting rational functions of the type chosen and 
with undetermined coefficients to the function approximated ; the coefficients are 
then determined to give a fit (in terms of the error function chosen to estimate 
the precision of the approximation) usually at the roots of a Chebyshev poly- 
nomial of suitable degree. The error functions are carefully plotted and the goal 
of the approximation is to ‘‘level’’ the error curves; that is, the coefficients and 
the form are adjusted, if necessary, to give an error function with the largest 
number of maxima and minima consistent with the approximation form chosen 
and with these maximal and minimal values all approximately equal in absolute 
value. The form chosen is frequently obtained by varying coefficients of TAYLOR’s 
series. 

Most readers unfamiliar with the work of the authors will be amazed by the 
precision obtained with a small number of coefficients. 
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The following approximations are polynomials not requiring divisions: 
arctan x (odd powers, three coefficients give 3D, eight coefficients give 8D), 
sin xx (odd powers, three coefficients give 4S, five coefficients give 8S), 10*, 
r(i + x), In (1 + x). Other approximations involve divisions. 

The early part of the book is devoted to an exposition of the techniques used 
by the authors in curve fitting. This is copiously illustrated—the ratio of space 
devoted to displayed illustration to that devoted to text appears to be about 5:1. 
The portion of this exposition which tells of the more or less subjective art of 
curve fitting as seen by these masters should be most valuable. Up to date no 
purely mechanical scheme has seemed to match their approximations, and their 
development of their method of curve fitting by example is most instructive. 
Earlier descriptions of the general properties of approximating functions are less 
impressive, for here summary statements in the form of theorems would appeal 
to the reviewer more than the narrative style used by the authors. However, 
bibliographic references at the ends of the chapters alleviate this objection. 

Use of these approximations and of these methods by the operators of high 
speed digital computers and punched card machinery is obvious. These machines 
have limited storage facilities for tables and incomparably greater computing 
facility for interpolation or for making the required approximation. Thus, for 
them these approximations, which are extreme in using complicated interpolation 
instead of listing more tabular values, are ideal. 

It would have been convenient to the operators of many machines to have had 
the coefficients listed in binary form also. The reviewer also hopes that future 
editions will contain a short list of iterative schemes commonly used in automatic 
computation for calculation of square roots, inverses, etc. The book as it stands 
is an aide which every coder or programmer of dan automatic machine will want ; 
with this slight additional information it would become indispensable. 

The type is easily read. Numbers of more than seven digits are divided into 
groups of four digits (the division might have been carried out for all numbers of 
more than four digits). Each error graph appears clearly on the same page as 
the approximation. 

cS ea 


1G. N. Warp, ‘‘The approximate external and internal flow past a quasi-cylindrical tube 
moving at supersonic speeds,”’ Quart. J. Mech. and Appl. Math., v. 1, 1948, p. 225-265. 


57[D, E, L].—W. F iucce, Four-place Tables of Transcendental Functions. 
McGraw-Hill Publishing Co., New York, 1954, 136 p., 23.5 K 15.5 cm. 
Price $5.00. 


Contents. sin x, cosx, tan x, and cot x, x = 0(.1)90°, cos x, sin x, tan x, 
cosh x, and sinh x, x = 0(.01)10.09 and tanhx, x = 0(.01)5.09. e* — 1 and 
1 — e-*,x = 0(.001).209, e* and e~*, x = 0(.01)10.09. In 10*, k = 1(1)10 and In x, 
x = 1(.01)10.09. Jo(x), Ji(x), Yo(x), Vilx), Io(x), I:(x), Ko(x), and K,(x), 
x = 0(.01)10.09. berx, beix, ber’ x, bei’ x, kerx, keix, ker’x and kei’ x, 
x = 0(.01)10.09. Elliptic integrals F(x, y) and E(x, y),x = 0(2)90° and y = 0(1)90°. 
erf x x = 0(.01)10.09, 1 — erf x x = 1(.01)3.09, Fresnel integrals C(x) and S(x) 
and exponential cos and sin integrals Ei(—<), Ei x, Cixand Sixx = 0(.01)10.09 
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and I'(x) x = 1(.001)2.09. All to 4D or 4S as appropriate. Also lists of formulas 
and a few notes concerning sources of more extensive tables and the use of these 
tables. In particular there are formulas for computing values beyond the range 
of the tables. 

The author does not indicate how the tables were computed or checked, but 
his bibliographic references indicate that they were ethically extracted from 
standard tables, to which he refers explicitly. He estimates errors at no more than 
0.6 in the decimal following the ones listed except for the Fresnel integrals, where 
his source material was admittedly weak. The Fresnel integrals are not in the 
notation of FLETCHER, MILLER, and ROSENHEAD (see 20.6 and 20.62 in their 


‘ gis 2 
Index) ; standard arguments are obtained by multiplying E by the square root 
T 


of the author’s arguments. (These standard arguments are used in both the Table 
of Fresnel Integrals by A. van WIJNGAARDEN and W. L. SCHEEN, Report R49 of 
the Mathematical Centre at Amsterdam, RMT 790, MTAC, vol. 4, p. 155-6, 
and the 1953 Tablitsy Integralov Frenelya of the Akad. Nauk SSSR, 40, MTAC, 
1955, v. 9, p. 75). A few values tested by the reviewer differed by at most 1 in the 
fourth decimal from values obtained by interpolation from the Akad. Nauk tables. 
The author mentions neither the Amsterdam nor the Akad. Nauk tables as pos- 
sible sources of values. 

The author seems to have put together a convenient set of tables for the use 
of engineers, the trigonometric tables combined with Bessel functions and other 
selected functions indicating fairly well the tables most likely to be reached for 
by a recently educated engineer. The tables will probably be popular, although 
this reviewer feels that in the long run an elementary set of tables (including some 
statistical tables) plus Jahnke and Emde would provide almost minimal require- 
ments for the expected users ; however there are some tables found in Fliigge’s book 
and not found in Jahnke and Emde—ker x and kei x, for example. 

hn Ba Be 


58[D ].—NBS Applied Mathematics Series No. 43, Tables of Sines and Cosines for 
Radian Arguments. U.S. Gov. Printing Office, Washington, 1955, xi + 278 p., 
26.7 X 20 cm. Price $3.00. 


Contains sin x and cos x, x = 0(.001)25.199, 8D; sin x and cos x, x = 0(1)100, 
8D; sin x and cos x, x = m X 10-?, m = 1(1)9, p = 1(1)5, 15D; sin x and cos x, 
= 0(.00001)0.00999, 12D ; conversion tables between radians, degrees, minutes, 
and seconds; and values of 4p(1 — p), p = 0(.001)1, exact. This is a new edition 
of table MT4, which was published in 1940 (RMT 81[D] MTAC, v. 1, 1943, 
p. 14-16. The present table has been extended .2 of a radian to finish a 4th whole 
cycle in the main table. A misprint in an argument in the early edition is noted, 
and this has been corrected in the present edition. Evidently no incorrect func- 
tional values have been noted, and no changes were reported in this edition. The 
last table in the book has been changed from a table of p(1 — p) to a table of 
3p(1 — p) over the same range of arguments p. 
The printing was designed for accuracy and economy rather than legibility. 
Most numbers are easily read. However, the publisher notes that the value for 
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cos 16.639 is hard to read in many volumes; this should be cos 16.639 = 
—0.5970 0260. 
CaF. 


59[F ].—I. M. Vrnocrapov, Elements of Number Theory. Translated from the 
fifth edition by Sau Kravetz. Dover Publications, Inc., New York, 1954, 
256 p. Paperbound edition $1.75. Clothbound edition $3.00. 

This volume contains two kinds of tables: 

(1) Tables of indices and powers of least positive primitive roots Modulo p 
for p < 100 (p. 220-225). These tables are equivalent to those in UsPENSKY and 
HEASLET [1 ], though the format is different. A description of their construction 
is given on page 112. A comparison with Uspensky and Heaslet reveals the fol- 
lowing misprints: 

p=71. For ind 56 = 43, read ind 56 = 42 
p = 97. For ind 11 = 6, read ind 11 = 86. 


(2) Table of primes <4000 and their smallest primitive roots (p. 226-227). 
Three errata are in this table: 


For p = 2763, read p = 2753. 


A comparison of this table against a SWAC calculation determined the following 
two errors: 

pb = 1459. Forg = 5, read g = 3 

pb = 3631. For g = 21, read g = 15. 


CHARLES A. NICOL. 
University of Texas 
Austin, Texas 


1J. V. Uspensky & M. A. HEASLET, Elementary Number Theory, New York and London, 
1939, p. 477-480. 


60[F].—ALBERT GLODEN, Table des Solutions de la Congruence x‘+1=0 
(mod p?) pour 4000 < » < 6000. 2 handwritten pages, 30 X 21.2 cm., de- 
posited in the UMT FILE. 


J’ai donné dans mon article “Sur quelques congruences d’ordre superieur,” 
[1], une méthode pour résoudre la congruence 


x‘ + 1 = 0 (mod p’). 


Dans cet article j’ai également publié les deux solutions minima pour p < 1000. 
La présente Table couvre I’intervalle 4000 < » < 6000. Elle donne les deux 
solutions x; et x2 qui sont < p?/2. On sait que les autres solutions < p* sont 
P — fe. P* — Xe. 
Cette Table sera étendue prochainement A toutes les valeurs de p < 10000. 
Je communiquerai mes résultats aux personnes qui en feront la demande. 
ALBERT GLODEN 


11, Rue Jean Jaurés 
Luxembourg 


1 Bulletin de la Société Royale des Sciences de Liége, 1950, p. 429-436. 
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61[F ].—ALBErt GLoDEN, Table des Solutions Minima de la Congruence x* + 1 =0 
(mod p*) pour 6000 < p < 10000. 2 handwritten pages, 26.7 X 20.3 cm., 
deposited in the UMT Fine. 


Celle Table est construite de fagon analogue 4 une Table précédente couvrant 
l’intervalle 4000 < » < 6000. Elle donne les deux solutions x, x. < p?/2. On sait 
que les autres solutions < p? sont p? — x, et p? — xo. 

Je communiquerai mes résultats aux personnes qui en feront la demande. 


ALBERT GLODEN 
11, Rue Jean Jaurés 
Luxembourg 


62(F ].—R. J. Porter, List of Irregular Determinants of Exponent 3n [80000 
< —D < 150000]. Typewritten manuscript of 72 pages, 20.4 X 12.3, on 
deposit in the UMT Fre. 


This list extends to D = — 150000 the two previous lists [UMT 155, MTAC, 
v. 7, p. 34] and [UMT 185, MTAC, v. 8, p. 96]. In the entire range up to 150000 
there are 5836 D’s of which 14 and 56 have exponents of irregularity 6 and 9, 
respectively. All others have exponents of 3. 
D. H. LEHMER 


University of California 
Berkeley, California 


63(G ].—Taxkayuxk1 Tamura, “Notes on finite semigroups and determination of 
semigroups of order 4,” Journal of Gakugei, Tokushima University, Japan, 
v. 5, 1954, p. 17-27. 


On pages 24 to 27 is a table of all 126 distinct semigroups of order 4. (For 
terminology see OLGA Taussky’s review in MTAC, v. 8, 1954, p. 231.) Each 
semigroup is represented by one multiplication table in the form 





AJaere 
gas o 
er tk QS 
ao oe & 











The semigroups are classified into six main categories : unipotent ; commutativity- 
indecomposable ; decomposable, 2-2 type ; decomposable, 3-1 type ; decomposable, 
2-1-1 type; commutative, non-uninpotent. 

The tables were computed by hand by M. YAMAMuRA and T. Tamura in 1953, 
and again in 1954, and were submitted for publication in August 1954. 

The reviewer used SWAC in April 1955 to verify that Tamura’s tables are 
equivalent (in the sense of isomorphism or anti-isomorphism) to the 126 semi- 
groups computed by the reviewer with SWAC in May and June 1954. (See the 
review cited above; the SWAC tables are to appear shortly in Proc. Amer. Math. 
Soc.) Thus Tamura’s table is undoubtedly complete and without misprints (except 
in the classifying labels 3.1-4 and 3.1-8). 








TI 
accur. 


Univer 
Los Ar 


64[H 
tic 
to 
Bi 

of ref 

the m 

rathe 

comp 
of lin 
those 
cially 








of 


ny 


or 








REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 127 


The creators of this table are to be congratulated on the completeness and 
accuracy with which they carried out a very long calculation. 


G. E. ForsyTHE 
University of California 
Los Angeles, California 


64[H ].—NBS Applied Mathematics Series No. 29, Simultaneous Linear Equa- 
tions and the Determination of Eigenvalues. U.S. Gov. Printing Office, Washing- 
ton, 1953, iv + 126 p., 26.7 X 20 cm. Price $1.50. 


This volume consists of a collection of 19 papers representing the majority 
of reports presented at a Symposium in Los Angeles on August 23-25, 1951. For 
the most part the papers are short and descriptive giving a general view of results 
rather than the full mathematical details which are adequately covered in ac- 
companying references. By bringing together various methods for the solution 
of linear equations and eigenvalue problems the book serves a useful purpose for 
those who are concerned with the explicit computation of such solutions, espe- 
cially by means of automatic computing machines. 

A list of the papers and a brief synopsis of each follows: 


1. Tentative classification of methods and bibliography on solving systems 
of linear equations, by GEORGE E. ForsyTHE. 

2. Simultaneous systems of equations, by A. M. OstrowskI. 

3. The geometry of some iterative methods of solving linear systems, by 
Aston S. HOUSEHOLDER. 

4. Solutions of linear systems of equations on a relay machine, by CARL-ErIk 
FROEBERG. 

5. Some special methods of relaxation technique, by EpUARD STIEFEL. 

6. Errors of matrix computations, by PauL S. Dwyer. 

7. Rapidly converging iterative methods for solving linear equations, by J. 
BARKLEY ROssER. 

8. Some problems in aerodynamics and structural engineering related to 
eigenvalues, by R. A. FRAZER. 

9. Inclusion theorems for eigenvalues, by H. WIELANDT. 

10. On general computation methods for eigenvalues and eigenfunctions, by 
GAETANO FICHERA. 

11. Variational methods for the approximation and exact computation of 
eigenvalues, by ALEXANDER WEINSTEIN. 

12. Determination of eigenvalues and eigenvectors of matrices, by MAGNUS 
R. HESTENES. 

13. New results in the perturbation theory of eigenvalue problems, by F. 
RELLICH. 

14. Bounds for characteristic roots of matrices, by ALFRED BRAUER. 

15. Matrix inversion and solution of simultaneous linear algebraic equations 
with the IBM 604 electronic calculating punch, by GzorcE W. Petrie, III. 

16. Experiments on the inversion of a 16 X 16 matrix, by Jon Topp. 

17. A method of computing eigenvalues and eigenvectors suggested by clas- 
sical results on symmetric matrices, by WALLACE GIVENS. 
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18. Computations relating to inverse matrices, by JacK SHERMAN. 
19. Results of recent experiments in the analysis of periods carried out in the 
Istituto Nazionale per le Applicazioni del Calcolo, by GAETANO FICHERA. 


Paper 1 classifies known methods of solving linear equations into six main 
groups, and lists a bibliography of some 500 titles with cross-reference to the 
methods. Paper 2 deals with nonlinear equations and essentially describes a 
method of obtaining an estimate of the error of an approximate solution. Paper 3 
shows how various standard iterative methods of solution can be regarded as 
particular aspects of a general geometric method of projections. Paper 4 presents 
experimental results on inverting matrices of order up to 20 on the Swedish relay 
machine BARK by the Gauss and Jordan elimination methods. Paper 5 describes 
a general relaxation method and a special form of it that permits the solution of 
a system of linear equations of order in m steps. Paper 6 discusses in a very 
general way the errors (inherent and computational) in solving linear equations 
or inverting matrices; it contains a bibliography of 30 papers. Paper 7 presents 
an orthogonalization method which generates a sequence of trial solutions whose 
residuals are mutually orthogonal relative to a positive definite matrix; the pro- 
cedure of paper 5 belongs to this method. Paper 8 describes ways of approximating 
certain continuous physical systems by discrete systems which are amenabie to 
computation. Paper 9 solves the following problem for a broad class of matrices: 
let x and y be fixed vectors, and consider the matrices A such that y = Ax; 
characterize those point sets in the complex plane which contain at least one eigen- 
value of every matrix A. Paper 10 presents first a method of finite least squares 
approximation for operators in Hilbert space and second a finite approximation 
method for differential equations based on the Cauchy-Lipschitz procedure. Paper 
11 derives inequalities for eigenvalues of a positive definite completely continuous 
operator on a subspace Q of a real Hilbert Space H by considering a sequence of 
intermediate problems corresponding to a chain of diminishing subspaces from 
H to Q. Paper 12 describes the following three methods for the solution of eigen- 
value problems: (i) power method, (ii) gradient method, (iii) an orthogonalization 
method closely related to that of paper 7. Paper 13 gives two results of E. Herz 
based on an essential inequality of Heinz. Paper 14 extends the standard ‘“‘circle”’ 
bounds for characteristic values to more general ‘‘oval’’ bounds. Paper 15 gives 
a program for carrying out the standard elimination method. Paper 16 summarizes 
tentatively the results of applying (i) the method of paper 15 on IBM equipment, 
(ii) the Monte Carlo method, and (iii) the iteration method X,,; = X,(2 — AX,) 
on the SEAC to a particular 16 X 16 matrix. Paper 17 uses the method of re- 
ducing a symmetric matrix to a triple diagonal form by a sequence of at most 
(1/2)(m — 1)(m — 2) planar rotations, m being the order of the matrix. Paper 18 
gives an algorithm for obtaining the inverse of B from the inverse of A when B 
differs from A in a column or in a row. Finally, paper 19 outlines an approximate 
method for finding the periods of an empirically given function on a finite interval. 

For the most part the papers employ mathematical methods which are ac- 
cessible to others than professional mathematicians and for this reason should 
interest a comparatively large audience to whom the problems treated are of 
importance. Paper 13, which employs the most advanced techniques of any of 
the papers, seems out of place in the series because of its purely theoretical con- 
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cern. For another reason, paper 8 strikes one as exceptional, for it expresses 
results mainly in physically descriptive terms leaving aside the mathematical 
problems of explicit calculation. But in the main, the papers bear directly on the 
apparent objective of the collection. 

Various of the papers contain numerical examples but from the point of view 
of one who wishes to compute solutions, especially with high speed machines, 
what is required is considerably more experimental evidence to distinguish the 
computational merits or faults of the proposed methods, particularly for systems 
of high order. Undoubtedly new information of this kind has been obtained since 
the time of the Symposium, and it would be of interest to have it available. 

WILiiAM KarusH 
University of Chicago 
Chicago, Illinois 


65[H ].—NBS Applied Mathematics Series No. 39, Contributions to the Solution 
of Systems of Linear Equations and the Determination of Eigenvalues. U. S. 
Gov. Printing Office, Washington, 1954, iii + 139 p., 26.7 X 20 cm. Price 
$2.00. 


This volume is a sequel to the one reviewed above and concerns itself with the 
same general problem. As the term, “Contributions,” in the title suggests, the 
separate papers in the collection are less closely connected than those in the 
preceding volume, but they are sufficiently united in general purpose to serve as 
a useful companion to the earlier collection. The papers vary in length from four 
pages to fifty-four pages and include such topics as practical methods for compu- 
tation, iterative methods of solution in Hilbert space, numerical experiments with 
matrices, condition of a matrix, and bounds for the rank and eigenvalues of a 
matrix. Descriptions of the individual papers now follow. 

The first paper, “‘Practical solution of linear equations and inversion of ma- 
trices,’’ (54 p.) by L. Fox describes various methods from the point of view of 
the computer using a desk calculator, and illustrates the methods by some 
twenty-six examples with matrices of order 6. It is a very convenient compilation 
which allows a ready comparison of various techniques. After a brief treatment of 
iterative methods of successive approximation the author turns to a group of 
elimination methods including Gauss elimination, JORDAN elimination, and 
AITKEN’s “below-the-line’”’ scheme. The next group of methods are compact 
elimination methods, among which are the procedures of DooLitTLE and Crovt. 
These elimination methods are related to the triangular decomposition A = LDU 
(L = lower diagonal unit matrix, U = upper diagonal unit matrix, D = diagonal 
matrix), and then attention is turned to schemes depending explicitly upon this 
decomposition. Among these are the procedures of CHOLESKY, DwyER-WauGH, 
BANACHIEWICZ, and the author. Following this, orthogonalization methods are 
described, and finally the complex matrix is singled out for special attention. 

The second paper, ‘“‘Punched-card experiments with accelerated gradient 
methods for linear equations,"’ (16 p.) by A. I. and G. E. ForsyTHE reports the 
results of numerical experiments on two matrices of order 6 designed to test some 
proposed devices to accelerate the method of steepest descent, which often con- 
verges too slowly for practical computation. The results were favorable to a pro- 
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cedure which inserts a certain accelerating step periodically in the iteration by 
steepest descent. 

The third paper, “Iterative methods of solving linear problems in Hilbert 
Space,” (34 p.) by R. M. Hayes, is a report on a doctoral dissertation dealing 
with the solution of Ax = 0 in Hilbert space for a certain class of linear operators 
A. After establishing general convergence theorems, the method is specialized 
for positive definite operators to the gradient method, the conjugate direction 
method, and the conjugate gradient method; iteration formulas for successive 
approximations are given and rate of convergence is studied. Finally, the theory 
is applied to ordinary differential equations, integral equations, and partial 
differential equations. The investigation is essentially theoretical and is closely 
related to the work of M. R. HESTENEs. It would be valuable to have numerical 
data available pertinent to the application discussed in the paper. 

The fourth and fifth papers are respectively, ‘‘Tables of inverses of finite 
segments of the Hilbert matrix,” (4 p.) by I. R. SavaGe and E. LuxKacs, and 
“The condition of the finite segments of the Hilbert matrix,” (8 p.) by Joun 
Topp. In the former paper a formula is derived for the elements of the inverse of 
the matrix H, = |\1/(¢ +7 — 1)||, 4,7 = 1,2, ---,m, and the inverse H,- is 
tabulated for m = 2, 3, ---, 10. In the latter paper various measures of the ‘‘con- 
dition” of a matrix are estimated for H,, corroborating that H, is ill-conditioned. 
Computations of H,,—' by a standard program were carried out on the SEAC and 
broke down at such low orders as nm = 6 and 7. 

The sixth paper, ‘“‘Lower bounds for the rank and location of the eigenvalues 
of a matrix,” (14 p.) by K. Fan and A. J. HoFFMAN deals with Problem 1 : deter- 
mination of lower bounds for the rank of an arbitrary matrix A = (a,;) of order n; 
and Problem 2: determination of p; such that every eyein d of A satisfies at 
least one of the inequalities |\ — ai:| < p:(i = 1,2, ---,m). The best-known 
criterion for the non-singularity of A is |a,;| > p> |ai;|. This criterion is used and 


generalized, as is the known solution p; = > lau! to Problem 2. Problems 1 and 
ji 


2 are closely related and this is clearly demonstrated in the paper. In the last 
part of the paper A is taken to be normal and inequalities involving k eigenvalues 
are established which generalize results of WIELANDT and of WALKER-WESTON 
for k = 

The last paper, ‘Inequalities for the eigenvalues of Hermitian matrices,”’ 
(9 p.) by K. Fan deals with the eigenvalues \; > 2 > --- > A, of an Hermitian 
matrix H. Sufficient conditions on a set of numbers d; are established for the 


inequalities 4; < d; (¢ = 1,2, ---,m) to hold; the conditions are expressed in 
terms of |a;;| and > |a,;|* (¢ = 1, 2, ---, m). In the second part of the paper the 
j>i 


d; are related to the eigenvalues k; > k2 > --- > k, of a block diagonal matrix 
obtained by replacing by 0 all elements of H except those belonging to principal 


submatrices; e.g., the inequality ya k< > Mi, (4 = 1,2,--+,m) is shown 


i=1 
to hold. 
WILLIAM KaruUSH 


University of Chicago 
Chicago, Illinois 
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66(H, P].—H. F. P. Purpay, Linear Equations in Applied Mechanics. Inter- 
science Publishers, Inc., New York, 1954, xiv + 240 p., 22.1 & 13.9 cm. 
Price $3.50. 


This is an elementary textbook in which the author determinedly faces the 
problem of obtaining numerical solutions to the equations he sets up. It is filled 
with examples and illustrations. Theorems and precise mathematical statements 
are almost completely absent. The book is essentially a set of recipes lucidly illus- 
trated by numerical calculations and applications to mechanics. The author states 
in his introduction that “this book is intended to explain the easier parts of a 
cross-section or slice of mathematical study; not so much a selection of isolated 
topics as a group of related studies. No attempt is made at logical rigour in the 
modern sense and the main object is to help the reader get acquainted with the 
main ideas, notations, and results by the easiest means, bearing in mind the final 
objective which is to obtain substantially correct numerical results to definite 
problems.” 

The author’s approach may be illustrated by quoting some advice he gives: 
if an equation has repeated roots make some trifling or negligible alterations to 
the constants thus separating the roots (like the Belfast constable who dragged 
a dead horse from Chichester Street to May Street to simplify his report). 

The bibliographical references are incomplete. The author states that a con- 
siderable body of literature not easily accessible has grown up around the subject 
of the simultaneous solution of linear algebraic equations. At least as regards the 
accessibility of the literature, this statement is misleading. 

The numerical examples seem to be correctly solved, and it is impressive and 
encouraging to see an author devote himself so directly to the task of arriving at 
numerical solutions. The references to mechanics and other engineering applica- 
tions are lucid, and they help the exposition. Although the methods used are 
sound, they are not always fully explained. Thus, in solving the heat equation by 
a difference equation approximation, the author takes his time increment to be 
small relative to the square of the size of his length increment, but he gives no 
real explanation of why he does so. In solving systems of algebraic linear equations 
by relaxation, the author gives no automatic program of relaxation and a careless 
reader might impute to him an implication that none is known. It seems likely 
that use of this book in a course in a university in the United States would 
require that the instructor furnish a great quantity of additional material— 
possibly in another text book. However, many students may be expected to benefit 
from the advice and the examples contained in this book—especially if they are 
made aware of the mathematical background of the methods of calculation used 
and of the existence of other methods. 

The book covers material in linear algebraic equations and matrices, ele- 
mentary loading problems on beams, elementary finite differences, numerical 
solution of ordinary differential equations (with no mention of any predictor- 
corrector process), scalar and tensor fields with examples from mechanics, par- 
tial differential equations, and integral equations. No nonlinear problems are 
mentioned. 

CG. B's 
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67[K, M, Z].—K. D. Tocuer, “The application of automatic computers to 
sampling experiments,’ Roy. Stat. Soc., Jn., v. 16, 1954, p. 39. 


This paper describes some of the problems that arise when automatic com- 
puters are used for conducting sampling experiments. The generation of random 
elements is discussed in detail and some methods of producing random variables 
with common distributions are described. The use of sampling methods to evaluate 
multivariate integrals is discussed. Finally a programme is devised to conduct a 
special form of a restricted random walk. 


J. WEGSTEIN 
National Bureau of Standards 
Washington 25, D. C. 


68[L, M].—G. Eason, B. Nose, & I. N. SNEDDON, “‘On certain integrals of 
Lipschitz-Hankel type involving products of Bessel functions,’’ Royal Soc. of 
London, Phil. Trans. Sect. A, Math. and Phys. Sciences, v. 247, 1955, p. 
529-551. 


Contains the following 11 tables all to 4D and for p, ¢ = 0(.2)2(1)3, where 
p = b/aand ¢ = c/a. 


Table 1. a? f ” Jo(at) Jo(bte-e'tdt. 
Table 2. a? f ” Jo(at) Ix (bt)e-e'tdt. 
Table 3. a? f ” Js(at)Jo(bt)e-ettdt. 
Table 4. a? f ” Jy(at) Is (bt)emettdt. 
Table 5. a f ” Jo(at) Jo (bt)e-tdt. 
Table6. a f ” J:(at) Jo(bt)e-*tdt. 
Table7. a f ” Jo(at) Js (bt)e-e'dt. 
Table 8. a f J, (at) J, (bt)e—*'dt. 
Table 9. f ” Jo(at) J, (bt)e-e'dt/t. 
Table 10. f ” Jy(at)Jo(bt)e-e'dt /t. 


Table 11. f J; (at) J, (bt)e-*'dt/t. 
0 
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Also, auxiliary tables to facilitate interpolation in the neighborhood of a 
singular point in each of these tables. 

The publication. includes an introduction indicating possible uses for the 
tables, development of formulas upon which the integrations are based, recurrence 
relations and brief instructions for use outside the range of the tables, the tables 
themselves, auxiliary tables, and a bibliography. 

CaF. 


69(L ].—K. I. McKenzie & M. Rotuman, Table of Bi'(+x) for x = 0(0.01)2 and 
Bi’(—x) for x = 0(0.01)10. i +12 mimeographed pages, 33 K 20.5 cm., 
deposited in the UMT FILE. 


8D tables with 3,,? of the derivative of the Airy integral of the second kind, 
Bi(x). The derivative was computed from a table of reduced derivatives of Bi(x). 
Tables of logio Bi(x) and Bi’(x)/Bi(x) are being prepared by the authors. 


M. RoTHMAN 
Northern Polytechnic 
Holloway, London N.7 
England 


70[P ].—CuHEsTER Snow, “‘Formulas for computing capacitance and inductance,” 
National Bureau of Standards Circular 544. U. S. Gov. Printing Office, Wash- 
ington, 1954, 37 figures, 69 p. Price $.40. 


This is a collection of explicit formulas for the computation of the capacitance 
between conductors having a great variety of geometrical configurations, the 
inductance of circuits of various shapes, and the electrodynamic forces acting 
between coils when carrying current. Formulas for skin effect and proximity 
effect in concentric cables and parallel wires are included. The formulas for the 
simpler configurations are expressed in terms of elementary functions; in the 
formulas for more complicated shapes Legendre functions, elliptic functions, and 
Bessel functions occur. 

Section 1 gives formulas for capacitance. Section 2 gives formulas for induc- 
tance and electromagnetic force, including a derivation of the general formulas 
from a vector potential. Section 3 treats frequency effects. Section 4 is a brief 
section of Legendre functions and in particular it shows the computation of 
Legendre functions of order n + 1/2 (m integer) by means of tables of elliptic 
integrals. Section 5 contains the derivation of some of the more difficult formulas, 
especially those derivations involving spheroidal, toroidal or other curvilinear 
coordinates. Section 6 is a bibliography of 29 numbers. 

A. E. 


71[Z].—VAcLav Hruska (editor), CesKosLovENSKA AKAD. Vip. Laborato 
matematickych stroji. Matematické Stroje (Mathematical Machines). 
Nakladatelstvi Cesk. Akad. Véd, Praha, 1953, 132 p. 25 cm. 


The booklet is the first issue, made in 1320 copies and entitled ““Mathematical 
machines,”’ of a new series entitled ‘Machines for processing information” 
(Stroje na zpracovdn{ informacf). The booklet is issued by the Laboratory of 
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Mathematical Machines of the Czechoslovak Academy of Sciences. Editor of the 
booklet is Dr. Vaclav HruSka. A principal man in the Laboratory seems to be 
Antonin Svoboda. Assistants credited with helping prepare the booklet are K. 
Bém, V. Cerny, H. Jerma¥, Z. Korvas, K. KriStoufek, J. Marek, J. Oblonsky, 
O. Pokornd, Z. Pokorny, J. Raichl, F. Svoboda, H. Snelerov4, M. St&rbovA, 
M. Valach, and V. Vy%in. Professor E. Cech is mentioned as a member of the 
Academy of Sciences interested in the work. 

The booklet is 132 pages long, with three page summaries in both Russian and 
English. There are seven chapters. 

“Chapter 1 describes the general character of a modern automatic computer” 
(from the summary). Basic concepts are defined, and symbols and flow charts are 
introduced, independent of any particular machine. The Czechoslovak automatic 
computer SAPO is mentioned next, with a simplified diagram but no mention of 
the type of memory, speed, or such interesting details. 

SAPO is a binary machine, with a memory of 1024 cells (numbered octally) 
of 32 bits each. The 32nd bit is always a parity check digit. The words are often 
divided into: one parity bit, 24 significant digits, sign of exponent, five-bit ex- 
ponent, and sign of number. Sometimes they are in a coded floating decimal form, 
with six significant decimals. 

There are 27 different three-address commands, including division and several 
logical commands. 

Chapters 3, 4, and 5 explain programming and coding, using three examples: 
computing cosines, geometrical ray tracing, and the Runge-Kutta method for 
solving differential equations. 

Chapters 6 and 7 describe (with pictures) the punched card machinery of the 
Laboratory, built by the state corporation ARITMA. They use a 90 column card, 
really 45 double columns. Holes are round, and everything looks like Remington- 
Rand equipment. There is a calculating punch seemingly comparable with the 
IBM 602; there does not seem to be any card-programming. An example deals 
with computing coordinates of points associated with the cross-sections of turbo- 
compressor blades. The errors of interpolation are discussed. 

One feels strongly the lack of any engineering specifications for the machines. 
I see no way to decide whether SAPO’s memory is a magnetic drum, or cathode- 
ray tube, or other. 

It is impossible to determine from the booklet whether SAPO exists only on 
paper, is being built, or is actually operating ; this is a usual situation with reports 
concerning computing machines. 

G. E. ForsyTHE 


University of California 
Los Angeles, California 


72(Z].—T. Pearcey, G. W. Hit, & R. D. Ryan, “The effect of interpretive 
techniques on functional design of computers,” Australian Journal of Physics, 
v. 7, 1954, p. 505. 


A number of computer codes were analyzed by counting the relative frequency 
of the various commands as they occur in complete programs and in interpretive 
function blocks. The results suggest that computers, suitably designed to use 
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interpretive methods, could store the standard function blocks in a fixed or semi- 
fixed storage system of rapid access. A relatively small additional amount of 
rapid-access erasable store would be required and large amounts of problem data 
would be held in slower-access backing store. A comparison of certain interpretive 
techniques with direct coding techniques shows interpretive techniques to be 
roughly twice as expensive in time and half as expensive in storage space as direct 
techniques. Based on this study, details of a functional machine design are sug- 
gested and the advantages of the resulting code system are discussed. 
J. WEGSTEIN 


National Bureau of Standards 
Washington, D. C. 


73[Z ].—R. J. K6nic, Introduction to FLAC Coding. Air Force Missile Test 
Center, Patrick Air Force Base, Florida, v + 114 p. 27 cm. 


This is primarily an introduction to FLAC programming and contains a certain 
amount of background information on modern high speed digital computers in 
addition to specific FLAC coding information. Among its sixteen operations 
FLAC has certain notable operations such as “‘decimal-binary conversion,” 
“‘binary-decimal conversion,’ ‘‘tape advance m words,” “‘tape reverse » words,” 
“tape hunt,” “‘tally,” “read out m words,” etc. 


” 46 


9” 46 


J. WEGSTEIN 
National Bureau of Standards 
Washington, D. C. 


74[Z ].—Los ALamos Screntiric LaBoratory, MANIAC. 


This is a 306 page programmer's manual for the MANIAC. It contains 
chapters on coding and flow diagrams, binary arithmetic, a simplified discussion 
of the various computer components, ‘‘descriptive’’ coding and subroutines, and 
operating procedures. It is liberally illustrated with examples and is another good 
example of a programming manual. 

J. WEGSTEIN 


National Bureau of Standards 
Washington, D. C. 


75[(Z ].—Manual of the SWAC Computing System, National Bureau of Standards 
Institute for Numerical Analysis, Los Angeles, California, viii + about 300 p. 
continuously revised. 29 cm. 


This work views the SWAC Computing System from a point midway between 
that of the engineer and the programmer. Thus it is neither a programming 
manual, nor an operation and maintenance manual, but includes most of the 
essential features of both. 

The manual is divided into five sections. These consist of (1) an introduction, 
(2) a detailed description of the ‘“hardware”’ of the machine itself, (3) a discussion 
of the order code and of standard coding procedures, and (4 and 5) a discussion 
of operating and testing procedures and of the service routines used. 

Most of the manual is written from an engineering standpoint. For example, 
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the description of the machine itself describes the physical location of each piece 
of equipment and its use in great detail. 

The sections on operating and coding provide enough information to enable 
one unfamiliar with the SWAC to make use of standard service techniques in 
programming. With this end in mind, the chief criticisms are: 


(1) The information concerned with preparation of data for input and the 
processing of output data is scattered throughout the manual and hence is 
difficult to follow. 

(2) It would be hard to use the manual as a reference, since there is no de- 
tailed index and much of the information is given in descriptive style. (However, 
this is unfortunately true of most manuals of this type). 


JoHN W. Cooper 
National Bureau of Standards 
Washington, D. C. 


76(Z ].—PROGRAMMING RESEARCH GrRouP, APPLIED SCIENCE Dtvision, INTER- 
NATIONAL BUSINESS MACHINES CORPORATION, “Specifications for the IBM 
mathematical FORmula TRANslating system, FORTRAN, International 
Business Machines Corp., New York, 1954. 43 p. 28 cm. 


This contains proposed specifications for a code which IBM plans to prepare 
which will permit the IBM 704 computer to directly accept mathematical prob- 
lems written in their usual mathematical form. Besides mathematical expressions, 
there are provisions for accepting input-output formulas (card reading formulas, 
tape reading formulas, print formulas, etc.), and formulas for all other information 
necessary in solving a mathematical problem. 


. J. WEGSTEIN 
National Bureau of Standards 
Washington, D. C. 


77(Z ].—ALBERT H. RUBENSTEIN, Editor, Coordination, Control, and Financing of 
Industrial Research, Proceedings of the Fifth Annual Conference on Industrial 
Research, June, 1954, with selected papers from the Fourth Conference, June, 
1953. King’s Crown Press, Columbia University, New York, 1955. 


Includes nontechnical chapters on Introduction to Computer Technology by 
C. B. Tompkins, Application of High Speed Computers to Research Problems 
by RicHarp F. CLippInceR, Automatic Data Reduction by G. TRUMAN HUNTER, 
and Operation of an Industrial Computing Facility by H. R. J. Groscn, extending 
through about 30 pages. 
c..R.. 7. 


TABLE ERRATA 


In this issue references have been made to errata in Review 55 and Review 59. 
244.—T. LarBLE, ‘“‘Héhenkarte des Fehlerintegrals,’’ Zeit. angew. Math. Physik, 
1951, p. 484-487. [MTAC, v. 6, 1952, p. 232.] 


T. Laible gives the first five complex zeros of the error function fo* e~“du 
which he obtained graphically to 3D. His 2nd and 3rd zeros do not deviate by 
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NOTES 137 


more than 2 units in the third decimal from the correct values, but his 3rd and 
4th zeros are given as 3.375 + 3.640i and 3.805 +- 4.0451 whereas they should be 
3.335 + 3.6461 and 3.769 + 4.0611, respectively. 

HERBERT E. SALZER 


Diamond Ordnance Fuze Laboratories 
Washington 25, D. C. 


NOTES 


International Conference on 
“Electronic Digital Computers and Information Processing” 


An international conference on “Elektronische Rechenmaschinen und In- 
formationsverarbeitung”’ will be held on October 25-27, 1955 at the Institut fiir 
Praktische Mathematik (IPM) at Technische Hochschule, Darmstadt/Germany. 
The conference is sponsored by GAMM (Gesellschaft fiir Angewandte Mathe- 
matik und Mechanik) and NTG-VDE (Nachrichtentechnische Gesellschaft im 
Verband Deutscher Elektrotechniker). 

Detailed information may be obtained from 


Prof. Dr. A. Walther 

Institut fiir Praktische Mathematik (IPM) 
Technische Hochschule 
Darmstadt/Germany. 


QUERIES AND REPLIES 


It is frequently possible to obtain information about tables and other aids to 
computation by publishing a query in this section; these should be addressed to 
the Chairman of the Editorial Committee. If information about tables is desired 
the question should contain an explicit definition of the function whose values 
are desired, the ranges of arguments of interest, and the precision required. The 
definition of the function might frequently be accomplished by a reference to the 
literature, and it should usually be accompanied by a decimal classification of the 
function according to the scheme of Lehmer’s Guide to Tables in The Theory of 
Numbers, Fletcher, Miller, and Rosenhead’s An Index to Mathematical Tables; 
Science Abstracts A, or another index as appropriate. 

It is probable that replies to some queries will be sent directly to the originator 
of the query; it is requested that the originator of each query relay replies to the 
Chairman of the Editorial Committee for publication. 


CORRIGENDA 


See the last line of the note in this issue, p. 118, ‘‘Remark on determination 
of characteristic roots by iteration,” by J. L. Brenner and G. W. Reitwiesner, 
which refers to MTAC, v. 8, p. 238. 

V. 9, p. 45, CORRIGENDA, the upper limit of integration should read x not x 
in both integrals. 
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CLASSIFICATION OF TABLES 


Arithmetical Tables. Mathematical Constants 
Powers 

Logarithms 

Circular Functions 

Hyperbolic and Exponential Functions 
Theory of Numbers 

Higher Algebra 

Numerical Solution of Equations 
Finite Differences. Interpolation 
Summation of Series 

Statistics 

Higher Mathematical Functions 
Integrals 

Interest and Investment 

Actuarial Science 

Engineering 

Astronomy 

Geodesy 

Physics, Geophysics, Crystallography 
Chemistry 

Navigation 

Aerodynamics, Hydrodynamics, Ballistics 


Calculating Machines and Mechanical Computation 
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